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In  this  work  a  numerical  model  simulating  the  aeroelastic  characteristics  of  a  flexible 

membrane  wing  is  presented.  The  use  of  the  Navier-Stokes  equations  as  the  fluid  dynamic 

model  in  the  present  membrane  wing  theory  is  a  substantial  departure  from  previous  work 

on  sail  aerodynamics  which  has,  almost  universally,  adopted  a  potential  based  description 

of  the  flow  field.  The  two-dimensional,  viscous  aeroelastic  problem  is  nondimensionalized 

and  a  set  of  six  basic  dimensionless  parameters  is  derived  which  govern  the  physical 

problem.  An  additional  parameter,  the  frequency  ratio,  is  proposed  as  a  meaningful 

parameter  for  characterizing  the  harmonically  driven  unsteady  aeroelastic  problem.  A 

numerical  procedure  is  developed  for  the  solution  of  the  coupled  aeroelastic  problem  and 

is  shown  to  yield  results  that  are  in  agreement  with  available  analytic  solutions  for  several 

appropriate  limiting  cases.  The  numerical  procedure  is  also  shown  to  satisfy  certain 

identities  as  dictated  by  the  fundamental  fluid  dynamic  conservation  laws.  The  role  of 

viscosity  in  membrane  wing  aerodynamics  is  investigated  using  the  numerical  model  for 

both  steady  and  unsteady  flows.  These  investigations  are  facilitated  by  distinguishing  three 

distinct  classes  of  problems  which  are  associated  with  limiting  cases  of  the  dimensionless 

parameter  set.  The  aerodynamic  characteristics  at  Reynolds  numbers  between  103  and  104 


are  shown  to  differ  substantially  from  those  predicted  by  a  potential  based  membrane  wing 
theory.  The  role  of  viscosity  is  shown  to  be  preeminent  in  the  harmonically  forced  unsteady 
flow  about  a  membrane  wing.  In  this  case  in  the  influence  of  viscosity  is  enhanced  since  the 
acceleration  and  deceleration  of  the  freestream  velocity  strongly  influences  the  separation 
and  reattachment  of  the  flow.  The  periodic  appearance  and  collapse  of  recirculation  zones, 
along  with  an  attendant  adjustment  in  the  membrane  configuration,  result  in  an  aeroelastic 
response  which  may  not  be  characterized  as  a  simple  harmonic  response  at  the  freestream 
forcing  frequency.  A  three-dimensional  potential  based  membrane  wing  element  is  also 
derived  and  is  shown  to  yield  results  that,  for  several  limiting  cases,  are  in  agreement  with 
available  analytic  solutions. 
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CHAPTER  1 
INTRODUCTION 

1.1  Overview  of  the  Sail  Analysis  Problem 

Marchaj,  in  his  second  book  on  the  science  of  sailing  (Marchaj  1979),  begins  the 
section  on  sail  design  with  the  following  comments. 


Despite  the  fact  that  mathematics,  computers  and  wind  tunnel  testing  are 
playing  an  increasing  part  in  the  designing  of  sails,  sailmaking  as  well  as  sail  tuning 
are  still  strongholds  of  art  based  on  a  hit-or-miss  technique  rather  than  on  science.  . 
. .  After  all,  unlike  the  aeroplane  wing,  which  can  be  regarded  as  a  rigid  structure  whose 
shape  is  unaffected  by  variation  in  incidence  and  speed,  sail  shape  is  a  function  of  both; 
in  which  the  shape  of  the  sail  affects  the  pressure  distribution  and  vice  versa,  in  a  rather 
unpredictable  manner,  (p.  500) 


These  comments  by  Marchaj  suggest  two  things  concerning  the  analysis  of  marine 
sails:  first,  the  behavior  and  performance  of  a  sail  is  governed  by  the  aeroelastic  interaction 
of  a  fluid  dynamic  field  and  a  deformable  surface ;  and  secondly,  as  a  result  of  this  interaction, 
as  well  as  the  presence  of  other  complexities,  analytic  and  computational  methods  which 
have  enjoyed  considerable  success  in  the  design  of  aircraft  have  not  as  yet  proven  useful  to 
sail  designers.  It  is  likely  that  the  second  observation  concerning  the  usefulness  of  analytic 
methods  to  sail  designers  will  not  always  hold  true. 

As  an  illustration  of  the  potential  usefulness  of  computational  methods  in  the  analysis 
of  sails  consider  the  membrane  wing  of  Fig.  1 ,  which  is  shown  operating  near  a  boundary 
such  as  the  surface  of  the  sea.  Since  the  sail  is  located  very  near  the  boundary  it  is  immersed 
in  the  shear  layer  adjacent  to  the  free  surface.  The  kinematic  inflow  conditions  to  the  sail  are 
consequently  nonuniform  and  may  also  be  unsteady  due  to  the  presence  of  wind  gusts. 
Furthermore,  since  the  sail  is  not  stationary  but  rather  has  some  forward  velocity,  the  relative 
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inflow  velocity  far  upstream  of  the  sail  will  vary  in  both  magnitude  and  direction  along  the 
sail  span  as  well  as  possibly  varying  in  time.  Computational  methods  provide  a  method  of 
simulating  such  a  complex  flow  environment  which  would  otherwise  be  nearly  impossible 
to  reproduce  experimentally. 


inflow  velocity 


membrane  wing 


free  surface  ▼ 


Figure  1.1.  Schematic  of  a  membrane  wing  of  finite  span  operating  near  a  free  surface. 


1.2  Review  of  the  Literature 


1.2.1  Two-Dimensional  Potential  Flow  Based  Models 

The  vast  majority  of  the  published  works  related  to  membrane  wing  aerodynamics 
have  made  several  simplifying  approximations  concerning  both  the  elastic  characteristics  of 
the  membrane  itself  as  well  as  the  nature  of  the  surrounding  flow  field.  Perhaps  the  most 
significant  of  these  simplifying  assumptions  is  that  the  fluid  dynamics  can  be  adequately 
described  by  a  potential  based  model  of  the  flow  field.  In  addition  to  the  almost  universally 
adopted  potential  flow  assumption,  the  additional  approximations  associated  with  thin 
airfoil  theory  -  small  camber  and  incidence  angle  -  are  also  often  made  and  the  membrane 
itself  is  generally  considered  to  be  inextensible.  A  comprehensive  review  of  the  work 


published  prior  to  1987  related  to  membrane  wing  aerodynamics  is  given  by  Newman 
(1987). 

The  analysis  of  membrane  wings  begins  with  the  historical  works  of  Voelz  (1950), 
Thwaites  (1961)  and  Nielsen  (1963).  These  works  considered  the  steady,  two  dimensional, 
irrotational  flow  over  an  inextensible  membrane  with  slack.  As  a  consequence  of  the 
inextensible  assumption  and  the  additional  assumption  of  small  camber  and  incidence  angle 
the  membrane  wing  boundary  value  problem  is  linearized  and  may  be  expressed  compactly 
in  nondimensional  integral  equation  form  as 

>      dHy/a) 

o 
where  y(x)  defines  the  membrane  profile  as  a  function  the  x  coordinate,  a  is  the  flow 
incidence  angle  and  Cj  is  the  tension  coefficient.  Equation  (1.1)  has  been  referred  to  as  the 
'Thwaites  sail  equation'  by  Chambers  (1966)  and  simply  as  the  'sail  equation'  by  Newman 
( 1 987)  and  Greenhalgh  et  al.  ( 1 984).  This  equation,  together  with  a  dimensionless  geometric 
parameter  e  which  specifies  the  excess  length  of  the  membrane,  completely  defines  the 
linearized  theory  of  an  inextensible  membrane  wing  in  an  steady,  inviscid  flow  field. 

Different  analytic  and  numerical  procedures  have  been  applied  to  the  sail  equation 
in  order  to  determine  the  membrane  shape,  aerodynamic  properties,  and  membrane  tension 
in  terms  of  the  angle  of  attack  and  excess  length.  In  particular,  Thwaites  (1961)  obtained 
eigensolutions  of  the  sail  equation  which  are  associated  with  a  wing  at  the  ideal  (singularity 
free)  angle  of  incidence.  Nielsen  (1963)  obtained  solutions  to  the  same  equation  using  a 
Fourier  series  approach  which  is  valid  for  wings  at  angles  of  incidence  other  than  the  ideal 
angle.  Other  more  recent  but  similar  works  are  those  by  Chambers  (1966),  Vanden-Broeck 
and  Keller  (1981),  Greenhalgh  et  al.  (1984)  and  Sugimoto  and  Sato  (1988). 

Various  extensions  of  the  linear  theory  have  appeared  in  the  literature  over  the  years. 
Vanden-Broeck  (1982)  as  well  as  Murai  and  Maruyama  (1980)  developed  nonlinear  theories 


valid  for  large  camber  and  incidence  angle.  The  effect  of  elasticity  has  been  included  in  the 
membrane  wing  theories  of  Jackson  (1983)  and  Sneyd  (1984)  and  the  effects  of  membrane 
porosity  have  been  investigated  by  Murata  and  Tanaka  (1989).  In  a  paper  by  de  Matteis  and 
de  Socoi  ( 1 986)  experimentally  determined  separation  points  were  used  to  modify  the  lifting 
potential  flow  problem  in  an  attempt  to  model  flow  separation  near  the  trailing  edge.  The 
effect  of  elasticity  and  porosity  were  also  considered  in  this  work. 

Comparisons  of  the  various  potential  flow  based  membrane  wing  theories  with 
experimental  data  have  been  reported  by  several  authors  including  Greenhalgh  et  al.  (1984), 
Sugimoto  and  Sato  (1988)  and  Newman  and  Low  (1984).  In  general,  there  has  been 
considerable  discrepancy  between  the  measurements  made  by  the  different  authors  which 
have  all  been  in  the  turbulent  flow  regime  at  Reynolds  numbers  between  105  and  106 
(Jackson  and  Fiddes  1993) .  As  a  result  of  the  discrepancies  in  the  reported  data  -  primarily 
due  to  differences  in  Reynolds  number  and  experimental  procedure  -  the  agreement  between 
the  potential  based  membrane  theories  and  the  data  has  been  mixed.  In  particular,  the 
measured  lift  is  in  fair  agreement  with  the  predicted  value  when  the  excess  length  ratio  is  less 
than  .01  and  the  angle  of  attack  of  attack  is  less  than  5°.  However,  even  for  this  restricted 
range  of  values  the  measured  tension  is  significantly  less  than  predicted  by  theory. 
Furthermore,  for  larger  excess  lengths  and  incidence  angles,  both  lift  and  tension  are  poorly 
predicted  by  the  theory.  Flow  visualization  studies  indicate  the  main  reason  for  the 
disagreement  is  the  existence  of  a  thick  boundary  layer  or  region  of  separated  flow  on  the 
membrane,  typically  near  the  trailing  edge.  It  has  been  noted  by  several  authors  that  the 
presence  of  viscous  effects  such  as  thick  boundary  layers  and  separation  regions  will 
overshadow  any  implications  associated  with  the  linearizing  approximations  made  by  thin 
wing  theory  (Newman  1987). 

1.2.2  Two-Dimensional  Viscous  Flow  Based  Models 

Although  the  importance  of  viscous  effects  on  membrane  wing  aerodynamics  has 
been  recognized  for  quite  some  time  (Nielsen    1963,  Newman  and  Low  1984),  very  few 


works  have  been  published  which  address  the  issue.  Jackson  and  Fiddes  (1993)  coupled  a 
boundary  layer  calculation  to  a  potential  flow  model  which  was  then  combined  with  an 
inextensible  membrane  model.  The  first  use  of  the  Navier-Stokes  equations  as  the  fluid 
dynamic  model  in  a  membrane  wing  theory  appears  to  be  the  work  of  Smith  and  Shyy  ( 1 993). 
In  this  work  an  incremental,  material  coordinate  formulation  of  the  elastic  membrane 
problem  was  coupled  with  a  body-fitted,  finite  volume  formulation  of  the  Navier-Stokes 
equations.  Results  from  the  viscous  flow  based  membrane  wing  model  were  compared  with 
a  potential  flow  based  membrane  wing  theory.  Although  the  results  were  in  qualitative 
agreement  with  the  flow  visualization  studies  given  by  Newman  and  Low  (1984)  a 
meaningful  quantitative  comparison  with  the  experimental  data  was  not  possible  due  to  the 
assumption  of  laminar  flow  in  the  model. 

1.2.3  Three-Dimensional  Potential  Flow  Based  Models 

The  development  of  a  three-dimensional  aeroelastic  membrane  wing  model 
introduces  several  complicating  factors  in  the  analysis  which  do  not  have  to  be  considered 
in  the  two-dimensional  analysis.  In  contrast  to  the  two-dimensional  case  where  the 
inextensible  approximation  is  often  valid,  the  effects  of  elasticity  are  preeminent  in 
three-dimensional  wings  where  spanwise  twist  and  trailing  edge  deflection  are  strongly 
influenced  by  the  elastic  properties  of  the  membrane.  Also,  even  within  the  context  of 
inviscid  flow,  the  tension  in  a  three-dimensional  membrane  is  no  longer  a  single  constant 
value  but  is  best  described  by  a  state  of  biaxial  tension  directed  along  lines  of  principal  stress 
(Jackson  1985).  Furthermore,  if  one  of  the  principal  tensions  vanishes  the  membrane  will 
compress  and  wrinkle  which  further  complicates  the  analysis.  Finally,  the  effect  of  viscosity, 
which  is  known  to  affect  strongly  the  aerodynamics  of  two-dimensional  membrane  wings 
under  many  conditions,  must  also  play  an  important  role  in  the  behavior  of  finite  span 
membrane  wings. 

Several  works  have  appeared  in  the  literature  which  have  considered  various 
idealizations     of    the     general     three-dimensional     membrane     wing     problem.     The 


two-dimensional  theory  described  in  the  previous  section  was  combined  in  a  stripwise 
fashion  with  lifting  line  theory  by  Nielsen  ( 1 963),  Sneyd  ( 1984)  and  Ormiston  ( 197 1 ).  Murai 
and  Maruyama  (1982)  developed  a  model  for  a  rectangular  planform  wing  assuming  the 
tension  in  the  chordwise  direction  to  be  zero.  Holla  et  al.  (1984)  also  investigated  a 
rectangular  membrane  wing  with  fixed  edges  assuming  a  state  of  uniform  biaxial  tension  to 
exist  in  the  membrane.  Kroo  (1986)  developed  a  numerical  method  for  arbitrary  planforms 
and  edge  conditions  which  was  based  on  a  strain  energy  minimization  formulation  of  the 
membrane  problem  and  a  vortex  lattice  treatment  of  the  potential  flow  problem.  Kroo  ( 1986) 
used  this  model  to  investigate  the  stability  of  hang  gliders.  A  summary  of  moderate  aspect 
ratio  membrane  wing  theories  is  given  by  Newman  (1987). 

The  numerical  model  developed  by  Jackson  (1985)  and  Jackson  and  Christie  (1987) 
is  probably  the  most  general  three-dimensional  model  presented  to  date  that  has  been  applied 
to  the  aeroelastic  analysis  of  marine  sails.  Jackson's  model  combines  a  vortex  lattice 
treatment  of  the  lifting,  potential  flow  problem  with  a  Lagrangian,  finite  element  based 
formulation  of  the  elastic  membrane  problem.  An  iterative  procedure  is  then  used  to  solve 
the  coupled  aeroelastic  boundary  value  problem.  Results  for  rectangular  and  triangular 
planform  sails  with  several  different  edge  boundary  conditions  are  presented  in  the  1985  and 
1987  papers.  The  issue  of  membrane  wrinkling  is  also  discussed  in  these  papers. 

1.3  Focus  of  the  Present  Work 

The  primary  focus  of  the  present  work  is  to  develop  a  computational  model  of  a 
two-dimensional  membrane  wing  which  fully  accounts  for  the  effect  of  viscosity  in  the  fluid. 
The  incorporation  of  a  viscous  fluid  dynamic  model  in  the  analysis  distinguishes  the  present 
work  from  previous  analytic  and  computational  work  on  sail  aerodynamics  in  which  the 
fluid  dynamics  were  characterized  by  a  potential  based  description.  Many  authors  have  noted 
that  flow  separation  is  often  observed  in  flow  visualization  studies  with  membrane  wings  and 
several  have  gone  on  to  say  that  the  existence  of  flow  separation  may  be  the  primary  cause 


for  the  disagreement  between  the  experimental  data  and  existing  membrane  wing  theories 
(Newman  and  Low  1984,  Newman  1987,  Sugimoto  and  Sato  1988). 

In  Chapter  2  the  governing  fluid  dynamic  conservation  laws  are  presented  along  with 
the  spatial  coordinate  based  equations  of  equilibrium  for  a  flexible,  linearly  elastic 
membrane.  The  boundary  conditions  on  fluid  pressure,  membrane  geometry  and  fluid 
velocity  which  couple  the  membrane  equilibrium  equations  to  the  fluid  dynamic 
conservation  laws  are  also  presented.  The  governing  equations  and  boundary  conditions  are 
then  nondimensionalized  to  determine  the  dimensionless  parameters  that  characterize  the 
unsteady,  viscous,  aeroelastic  membrane  wing  problem. 

In  Chapter  3  the  discrete  form  of  the  governing  equations  is  given.  The 
Navier-Stokes  equations  are  first  transformed  from  Cartesian  to  general  curvilinear 
coordinates  and  a  pressure  based  control  volume  formulation  is  adopted  as  the  basic  point 
of  departure  for  the  solution  algorithm.  Special  attention  is  given  during  the  development 
of  the  numerical  method  to  the  consistent  transformation  of  the  components  of  the  fluid 
velocity  vector  from  Cartesian  to  curvilinear  coordinates.  Special  care  is  also  given  to  the 
numerical  treatment  of  the  time  derivative  term  in  the  transformed  continuity  equation  to 
insure  a  consistent,  source  free  implementation  of  mass  conservation.  The  membrane 
equilibrium  equations  are  implemented  in  a  straightforward  way  using  conventional  finite 
difference  approximations. 

The  numerical  method  developed  in  Chapter  3  is  applied  in  Chapter  4  to  several 
elementary  test  problems  for  which  analytic  solutions  are  known.  The  test  cases  presented 
for  the  Navier-Stokes  algorithm  are  chosen  to  illustrate  the  need  for  special  care  in  certain 
aspects  of  the  implementation.  The  test  case  given  for  the  membrane  equilibrium  algorithm 
is  included  simply  as  a  check  on  the  accuracy  of  the  computation. 

In  Chapter  5  the  computational  procedure  is  applied  to  three  separate  classes  of 
aeroelastostatic  membrane  wing  problems.  These  three  classes  of  problems  arise  naturally 
from  the  nondimensionalization  of  the  governing  equations  and  are  referred  to  as  the 
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constant  tension  case,  the  elastic  case,  and  the  inextensible  case.  Results  from  the  viscous 
flow  based  model  are  compared  with  results  using  a  potential  flow  model  of  the  fluid 
dynamics.  A  very  limited  comparison  of  the  computational  results  with  selected 
experimental  data  is  also  made  in  this  section.  A  comprehensive  comparison  with 
experimental  data  is  precluded  in  the  present  work  by  the  assumption  of  laminar  flow  since 
the  experimental  data  were  obtained  at  Reynolds  numbers  in  the  turbulent  regime. 

In  Chapter  6  the  numerical  method  is  applied  to  membrane  wings  subjected  to  a 
harmonically  varying  freestream  velocity.  The  three  classes  of  wings  are  investigated  and 
the  time  dependent  aerodynamic  properties  of  the  wing  are  compared  with  the  steady  state 
values  under  comparable  conditions. 

In  Chapter  7  membrane  wings  of  finite  span  are  considered  and  an  incremental, 
velocity  stepping,  finite  element  based  formulation  is  presented.  The  fluid  dynamics  are 
approximated  by  a  potential  based  description  of  the  flow,  and  a  material  coordinate  based 
description  of  the  membrane  kinematics  is  adopted.  A  three-dimensional  aeroelastostatic 
membrane  wing  element  is  derived  and  tested  on  several  limiting  cases  of  wing  geometry, 
pressure  distribution  and  material  stiffness  before  the  algorithm  is  applied  to  an  elliptical 
planform,  edge  constrained  membrane  wing  of  moderate  aspect  ratio. 

Finally,  in  Chapter  8,  the  present  work  is  summarized  and  several  observations  are 
made  concerning  the  role  of  viscosity  in  the  steady  and  unsteady  flow  about  a  membrane 
wing.  Suggestions  for  future  work  are  made  in  this  chapter  with  particular  emphasis  on 
incorporating  the  effect  of  turbulence  in  the  membrane  wing  model.  The  importance  of 
comparing  the  computational  predictions  with  experimental  data  is  also  discussed  in  this 
chapter. 


CHAPTER  2 
GOVERNING  EQUATIONS 

2. 1  Membrane  Equilibrium 

In  this  section  the  general  equilibrium  equations  are  presented  for  a  two-dimensional 
elastic  membrane  subjected  to  both  normal  and  shearing  stresses.  The  membrane  is  assumed 
to  be  massless  and  the  equilibrium  conditions  are  stated  in  terms  of  the  instantaneous  spatial 
Cartesian  coordinates  jc,  and  the  body-fitted  curvilinear  coordinates  £j.  The  basic 
formulation  is  essentially  identical  to  many  previously  published  works  such  as  de  Matteis 
and  de  Socoi  (1986)  and  Sneyd  (1984).  Index  notation  with  the  summation  convention  is 
used  throughout  the  following  chapters. 


membrane  free  body 


Figure  2.1.  End  constrained  elastic  membrane 


Figure  2. 1  illustrates  an  elastic  membrane  restrained  at  the  leading  and  trailing  edges 
subjected  to  both  fluid  dynamic  pressure  and  shear  stress,  p  and  x  ,  respectively.  Imposing 
equilibrium  in  the  normal  and  tangential  directions  requires 
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"V1+<&r_-/4rt  (2.1, 


dx}\         dxj  \T 


dT 

Wx    =   "  '  (2'2) 

where  7  is  the  membrane  tension.  The  net  pressure  and  shear  stress  acting  on  a  segment  of 
the  membrane  are  given  respectively  by 

Ap  =  p~  -  p+  (2.3) 

t  =  r~  +  r+  (2.4) 

where  the  superscript  indicates  the  value  at  the  upper  and  lower  surface  of  the  membrane  as 
shown  in  the  figure.  If  the  membrane  material  is  assumed  to  be  linearly  elastic  the  nominal 
membrane  tension  T  may  be  written  in  terms  of  the  nominal  membrane  strain  5  as 

T  =  (°S  +  Ed)h  (2.5) 

where  °S  is  the  membrane  prestress,  E  is  the  elastic  modulus,  and  h  is  the  membrane 
thickness.  The  nominal  membrane  strain  is  given  by 

d  =  ^-^  (2.6) 

where  Lo  is  the  prestrained  length  of  the  membrane  and  L  is  the  length  of  the  membrane  after 
deformation  which  may  be  expressed  in  terms  of  the  spatial  coordinates  x,  as 


L  = 


(2.7) 


where  c  is  the  chord  length. 

At  the  leading  and  trailing  edges  of  the  membrane  the  following  boundary 
conditions  are  imposed  on  Eq.  (2.1) 

x2  =  0    at    x,  =  0,c  (2.8) 
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2.2  Fluid  Dynamic  Conservation  Laws 

The  governing  conservation  laws  for  unsteady,  laminar,  incompressible  flow  are  the 
Navier-Stokes  equations  which  may  be  written  in  two-dimensional  Cartesian  coordinates  as 

|(  PV,  )  +  ^  (  ey,  )    =  JL  (  ^  )   -  JL  (2.9, 

|(  ev2  )  +  A  (  9yf2  )   -  JL  (  ^y  )  -  j|  (2.10) 

t  +  ^(ev;)=o  (2.ii) 

where  v;  are  the  Cartesian  components  of  the  fluid  velocity  vector,  q  is  the  fluid  density,  u, 
is  the  fluid  viscosity,  t  is  time,  and  p  is  the  fluid  pressure. 

When  new  independent  variables  %,\  and  ^2  ^Q  introduced  Eqs.  (2.9)  through  (2.11) 
change  according  to  the  general  transformation  |i=  ^1  (x^  X2,  t)  and  %2  -  I2  (xi,  x2>  0- 
Equations  (2.9)  through  (2.11)  can  be  rewritten  as  follows  where  the  subscript  i,j  indicates 
the  partial  derivative  of  the  i  Cartesian  component  of  velocity  or  position  with  respect  to  the 
general  curvilinear  coordinate  j. 

ft(  JQV,    )  +  |(^/,)     =   ^(  *f <  llV,j   "  ftV,;  ))  + 

a|"(  7<  ft"i;  "  <72"u  ))  "  ^-(^2,2  P)  +  ^-(^2.1  P)  (2.12) 

£(  /pv2  )  +  ■^■(QVjVj    =  ■£-(  fj(  qlVu  -  q2v22  ))  + 

a§-<  j(  ^2.2  -  92v2,i  ))  +  a|-(*u  p)  -  a|-(*u  p>  (2-i3) 


&(*>  +  4;(<^)  =0 


(2.14) 


where  the  contravariant  velocity  components  are  given  by 

V,    =   (v,  -J,)  *22   -  (v2  -*2)*i,2  (2.15) 

V2    =   (v2  -  *2)  x,  ,    -   (vj  -  x{)  x2A  (2.16) 
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and 


1\    =   (-ri.2)2   +   (*2.2)2  (2-17) 


?2    =   *1,1*1,2    +   *2,1*2,2  (2.18) 

<?3    =   (*i.i)     +   (*2,i)  (2.19) 

and  the  Jacobian  of  the  transformation  is  defined  as 

■*    =  ^l.l^^   ~~  X\2X2,\  (2.20) 

The  following  kinematic  boundary  conditions  are  imposed  on  the  Cartesian 
components  of  the  fluid  velocity  vector  at  the  membrane  surface 

v«  =  *i  (2.21) 

and,  in  the  present  work,  the  fluid  velocity  far  upstream  of  the  membrane  is  chosen  to  be  a 
harmonic  function  of  time  given  by 

v,  =  Voocosa(l  +/?sina>/  )  (2.22) 

v2  =  Va,sina(l  +  ft  sin  cot  )  (2.23) 

where  oo  is  the  circular  frequency  of  the  freestream  velocity  and  P  is  a  parameter  defining 
the  magnitude  of  the  oscillation. 

2.3  Nondimensionalization  of  the  Governing  Equations 

The  aeroelastic  boundary  value  problem  can  be  written  in  nondimensional  form  after 
introducing  the  following  dimensionless  variables 

vi 

vl  =  VZ  (2-24) 


h  =  £  (2-25) 


t  =  too 


(2.26) 


*1 

*1  =  T  (2-27) 
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x2 


x2  =  -f  (2.28) 


An 
Ap  =  -^-  (2.29) 

t  =  ^-  (2.30) 

°Sh 


or 


where  either  Eq.  (2.30)  or  Eq.  (2.31)  is  used  to  nondimensionalize  the  membrane  tension 
depending  on  whether  the  tension  is  dominated  by  pretension  or  by  elastic  strain. 
Substituting  these  variables  into  Eq.  (2.1)  leads  to  the  following  dimensionless  equilibrium 
equation 

when  membrane  tension  is  dominated  by  elastic  strain,  with  77/  defined  to  be 

*>  -  (£i 

or 

when  membrane  tension  is  dominated  by  pretension,  with  TIi  defined  to  be 

U2    =   (-?&)  (2.35) 

\Qvl,c) 

The  use  of  the  cube  root  in  the  definition  of  77;  in  Eq.  (2.33)  is  suggested  by  the  exact  solution 
of  Eq.  (2. 1 )  given  by  Seide  ( 1 977) . 

The  boundary  conditions  for  the  membrane  equilibrium  equation  may  also  be  written 
in  nondimensional  form  as 
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x2  =  o    at    .r,  =  0,  1  (2.36) 

where  the  hat  has  been  dropped  for  convenience  from  the  dimensionless  variables  in  Eqs. 
(2.32),  (2.34),  and  (2.36). 

The  physical  significance  of  the  aeroelastic  parameters  77/  and  IJ2  is  that  the 
nondimensional  deformation  of  an  initially  flat,  elastic  membrane  is  inversely  proportional 
to  77;  in  the  absence  of  pretension.  Alternatively,  the  deformation  of  a  membrane  is  inversely 
proportional  to  772  in  the  presence  of  large  initial  pretension.  Consequently,  the  steady  state, 
inviscid,  aeroelastic  response  of  an  initially  flat  membrane  wing  at  a  specified  angle  of  attack 
is  controlled  exclusively  by  77/  in  the  limit  of  vanishing  pretension  and  exclusively  by  77? 
in  the  limit  of  vanishing  material  stiffness. 

Substituting  the  same  dimensionless  variables  into  Eqs.  (2.9)  through  (2.11)  leads 
to  the  following  nondimensional  form  of  the  incompressible  Navier-Stokes  equations 

5'l<".  >  +  s;<v.»  =  (i)  ^  <"'■>»- I;         cud 

jgr  (  Vj  )    =  0  (2.39) 

with  boundary  conditions  at  the  membrane  surface 

v,  =  St  xt  (2.40) 

and  far  upstream  from  the  membrane 

V!  =  cosa(l  +£sinr  )  (2.41) 

v2  =  sina(l  +  /Ssinf  )  (2.42) 

where  again  the  hat  has  been  dropped  for  convenience  from  the  dimensionless  variables.  The 
dimensionless  parameters  appearing  in  Eqs.  (2.37)  and  (2.38),  the  Reynolds  number  and  the 
Strouhal  number  are  defined  as 
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_  VXQC 

Rn   =  —fT~  (2-43) 

St  =  !§£■  (2.44) 

If  the  membrane  is  not  initially  taut  the  geometry  of  the  wing  may  be  characterized 

by  an  additional  dimensionless  quantity,  the  excess  length  e,  defined  as  follows. 

.  Lq-c  (2.45) 

t  c 

In  order  to  relate  the  membrane  excess  length  to  a  more  familiar  airfoil  property, 

Greenhalgh  et  al.  (1984)  offered  the  following  approximate  formula,  based  on  a  parabolic 

arc  airfoil,  which  relates  e  to  the  wing  camber  by 

where  the  left  hand  side  of  Eq.  (2.46)  is  taken  to  be  the  value  of  the  X2  coordinate  at  the 
midchord  point.  The  set  of  dimensionless  parameters  given  above  -77;,  TT^,  St,  Rn,  e,  and 
a  -  completely  characterize  the  physical  problem  considered  in  the  present  work. 

In  addition  to  the  basic  dimensionless  parameter  set  other  physically  significant 
dimensionless  parameters  may  be  derived.  In  particular,  a  frequency  ratio  parameter,  Q,  may 
be  introduced  by  drawing  an  analogy  between  a  one  degree  of  freedom  spring/mass  system 
and  the  transverse  motion  of  the  membrane  in  response  to  a  harmonically  varying  freestream 
velocity.  If  the  ratio  of  the  system  forcing  frequency  to  the  system  natural  frequency  is 
defined  as 

Q=§-n  (2.47) 

the  following  parameters  may  be  defined  and  substituted  for  the  Strouhal  number  in  the  basic 
parameter  set. 

l 
Q2  =  jp-l  (2-49) 
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Again,  as  in  the  case  of  steady  flow,  Qi  is  the  appropriate  dimensionless  frequency  when  the 
membrane  is  pretensioned,  whereas  Q\  is  the  appropriate  dimensionless  frequency  when  the 
membrane  develops  tension  elastically. 

Finally,  the  aerodynamic  lift,  drag  and  moment  about  the  quarter-chord  as  well  as  the 
membrane  tension  and  fluid  dynamic  pressure  are  nondimensionalized  in  the  customary  way 
as  stated  below. 

(2.50) 
(2.51) 
(2.52) 
(2.53) 

(2.54) 

In  Eqs.  (2.50)  through  (2.52),  the  lift,  drag,  and  quarter-chord  moment  are  obtained 
by  integration  of  the  aerodynamic  pressure  and  shear  stress  over  the  membrane  chord  and 
resolving  the  net  force  vector  into  components  parallel  and  perpendicular  to  the  freestream 
velocity  vector. 


^L  ~ 

\qv\>c 

cD  = 

D 

{qv2*c 

c     — 

Mc/4 

*-M  ~ 

{qv^c2 

CT  = 

T 

jQv^c 

cp  = 

P  -  Pec 

lnv2 

CHAPTER  3 
NUMERICAL  METHOD 

3.1   Membrane  Equilibrium 

A  discrete  form  of  the  elastic  membrane  boundary  value  problem  can  be  obtained 
at  a  finite  number  of  points  on  the  fixed  interval  [0,c]  by  replacing  the  derivatives  in  Eq.  (2. 1) 
and  the  integral  in  Eq.  (2.7)  with  appropriate  finite  difference  and  finite  sum  approximations. 
Applying  central  difference  approximations  to  Eq.  (2.1)  leads  to  a  three  point  difference 
kernel  centered  around  point  P  with  neighboring  points  E  and  W  as  shown  in  Fig.  2. 1 .  At 
each  point  P  an  equation  of  the  general  form 

^P*2r  =  ^EX2E  +  Ay/X2w  +  Sxj  (3-1) 

is  then  obtained  where  the  As  are  coefficients  associated  with  the  finite  difference 
approximation,  and  5^2  is  a  source  term  containing  all  terms  in  Eq.  (2.1)  that  cannot  be 
expressed  as  a  linear  combination  of  the  xj  coordinates.  Consequently,  all  of  the  nonlinearity 
in  the  boundary  value  problem  is  contained  in  the  source  term.  The  resulting  set  of  finite 
difference  equations  may  then  be  solved  using  a  line  iteration  method  with  underrelaxation. 
As  a  measure  of  the  degree  to  which  the  discrete  form  of  Eq.  (2.1)  has  been  satisfied  the 
residual  of  the  set  of  membrane  equilibrium  equations  is  defined  as  follows. 


RXl  =    }_i\~  ApX2n  +  A.ex2e  +  Awx2w  +  SX2\  c 


(3.2) 


all  P 


3.2  Fluid  Dynamic  Conservation  Laws 

A  pressure  based  numerical  procedure  originally  proposed  by  Patankar  and  Spalding 
( 1 972)  for  Cartesian  coordinates  was  chosen  for  computing  the  laminar  incompressible  flow 
surrounding  a  two-dimensional  wing  of  vanishing  thickness  in  an  unbounded  domain.  The 
details  of  the  basic  pressure  correction  algorithm  are  given  in  Patankar  (1980)  with  the 
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extension  of  the  procedure  to  general  curvilinear  coordinates  given  in  Shyy  (1994).  The 
present  implementation  of  the  algorithm  follows  the  work  of  Braaten  and  Shyy  ( 1986)  with 
the  statement  of  the  conservation  laws  in  terms  of  a  time  dependent  grid  taken  from 
Thompson  et  al.  (1985). 


Figure  3.1.  Staggered  grid  arrangement  in  the  physical  domain 


A  staggered  grid  arrangement  is  adopted  for  the  discretization  of  Eqs.  (2.9)  through 
(2.11)  as  shown  in  Fig.  3.1.  Specifically,  a  nonorthogonal  body-fitted  grid  system  is 
generated  numerically  at  the  positions  marked  by  triangles.  Both  the  Cartesian  and 
contravariant  velocity  components,  vj  and  Vj  as  well  as  the  Cartesian  components  of  the  grid 
velocity  vector,  xt ,  are  located  at  the  midpoint  of  the  control  volume  faces  as  shown  in  the 
figure.  The  discrete  value  of  pressure  is  located  at  the  arithmetic  center  of  the  four  adjacent 
grid  points.  The  finite  volume  form  of  the  conservation  laws  may  be  obtained  by  integrating 
Eqs.  (2.9)  through  (2.11)  using  a  first  order,  fully  implicit  time  integration  scheme  over 
appropriately  staggered  control  volumes  leading  to  a  five  point  difference  kernel  centered 
around  the  point  P  with  neighbor  points  N,  S,  E,  and  W.  By  arbitrarily  taking  the  grid  spacing 
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in  the  transformed  domain  to  be  unity  and  the  time  step  to  be  A  t,  the  momentum  and  mass 
conservation  laws  take  on  the  following  discrete  control  volume  form 

(Jqvx-  J°qv0  ) 

-. — + 

At 

u  M 

[QVlV1  -  j(q\Vxl  -  q2vX2)  +  x22p)e  -  (qVxvx  -  -j(g,v,  ,  -  q2vx2)  +  x22p)w  + 

QV2V\   -  j(~  <?2V1,1   +  <?3V1,2)  ~  x2,\P">n  ~  (QV2v\   ~  /("  <?2V1.1   +  ^3V1.2>  ~  x2,\P)s 

=  0  (3-3) 

(   JQV2  ~  J°QV°.   ) 

A ~   + 

At 

'-QVlv2  ~  j(4\v2,\   ~  <?2V2,2)  +  x\3P)e  ~  ^Vlv2  ~  y<^lv2,l   _  ^2V2,2)  +  *l,2P)*f  + 
£V2V2  _  y(-  <?2V2,1  +  ^3V2,2)  -  *l,l/>)n  "  ^V2V2  ~  y(-  ?2V2.1   +  <?3V2,2)  ~  x\jP>» 
=  0  (3.4) 


(7g  Jr/)g)  +  (ev,),  -  (pv,)  „  +  (ev2)B  -  (ev2),  =  o  (3  5) 


where  the  superscript  refers  to  the  previous  time  level.  The  Cartesian  components  of  the  grid 
velocity  vector  appearing  in  Eqs.  (2. 1 5),  (2. 1 6)  and  (2.21)  are  approximated  by  the  first  order 
backward  time  difference  given  below. 


X:  -  x° 


X;    = 


At  (3.6) 

In  order  to  clarify  some  of  the  details  of  the  numerical  implementation  consider  the 
control  volume  used  to  derive  the  discrete  form  of  the  v/  momentum  equation  as  shown  in 
Fig.  3.2. 
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Figure  3.2.  Control  volume  used  for  the  v/  momentum  equation 


If  the  convective  momentum  flux  is  approximated  using  a  first  order  upwind  scheme 
and  the  diffusive  flux  is  approximated  by  a  central  difference  expression  the  discrete  form 
of  the  vj  momentum  equation  may  be  written  as 


APV\P  =  AEV\E  +  AwV\w  +  anv\„  +  Asvis  +  Sv, 


(3.7) 


where,  after  explicitly  imposing  mass  conservation,  the  coefficients  in  Eq.  (3.7)  are  given 
by 


AE  = 


Aw  — 


AN  ~ 


AS  = 


-  qvxI 

,  o 

qv{\w  , 

0 

-  QV2\n 

,  o 

gv2\s , 

0 

+  J1ll 


1./*    I 


+  7^3'" 


+  7^ 


AP  ~  AE  "*"  AW  +  AN  **"  AS  ~*~    jf 


(3.8) 

(3.9) 

(3.10) 

(3.11) 
(3.12) 


with  the  source  term  given  by 
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U  I  U  I  fi  I  jU  I 


7°pv0 


+  x22p\w  -  x22p\e  +  x2Ap\n  -  x2lp\s  +  —^ 


(3.13) 


where  the  various  terms  may  be  evaluated  at  the  control  volume  faces  using  appropriate 
interpolations.  In  the  present  work  the  second  order  upwind  scheme  of  Shyy,  Thakur  and 
Wright  (1992)  is  adopted.  As  a  result  of  the  higher  order  interpolation  scheme  used  for  the 
convection  terms,  additional  source  terms  will  appear  in  Eq.  (3.13).  The  explicit  form  of 
these  additional  terms  may  be  found  in  that  reference. 
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Figure  3.3.  Control  volume  used  for  the  V2  momentum  equation 

Similarly,  the  discrete  form  of  the  V2  momentum  equation  may  be  derived  by 
considering  the  control  volume  shown  in  Figure  3.3.  The  discrete  form  of  the  equation  may 
be  written  as 


ApVj    =  AFv2    +  Awv2    +  ANv2    +  Asv2   +  Sv 


^Ev2 


Sv2< 


(3.14) 


with  the  coefficients  given  by 
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AF  = 


-  gvAe  ,  o 


+  5*il. 


(3.15) 


Ajy    - 


ev,L  .  oj  +  ^,l 


t* 


(3.16) 


A^  — 


-^2I»   '    Oj+J^'n 


/" 


(3.17) 


A5  =  [  QV2\S  , 


+  1^ 


Ap    —    r\.  p    T   /i  yy     i     A^   ~r    a^    T 


and  the  source  term  given  by 


U  .  pL  \  M  I  f*  I 

5V2  =    -  jq2V2j}e  +  7?2v2,2'>v  ~  y^.l'n  +  y^^l'j 

yOpV0p 

+  X,  2ple  -  Jflt2plw  -  *upl„  +  Xup\s  +      At  P 


(3.18) 

(3.19) 


(3.20) 


Similarly,  the  discrete  form  of  the  pressure  correction  equation  may  be  derived  by 
considering  the  control  volume  shown  in  Fig.  3.4. 


P  w 


P    E 


Figure  3.4.  Control  volume  used  for  the/?'  equation 
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The  discrete  form  of  the  pressure  correction  equation  may  be  written  as 

AP  p'P  =  AE  p'E  +  Aw  p'w  +  AN  p'N  +  As  p's  +  Sp,  (3.21) 


with  the  coefficients  given  by 


AE  =  Q 


Aw  =  Q 


2,2 


I      " 


+ 


'1,2 


'x2  X2    ^ 

At         Ai 


(3.22) 


(3.23) 


AN  =  Q 


X" 


+ 


2,1 


(3.24) 


As  =  Q 


x2        x2  1 

Al,l         A2,l 


(3.25) 


AP  =  AE  +  Aw  +  AN  +  As 


(3.26) 


and  the  source  term  given  by 


s,  =  q(K  ~  K  +  vi  -  vi)  +  ^7^ 


(3.27) 


The  sequence  of  events  in  the  solution  procedure  is  as  follows.  The  momentum 
equations  are  first  solved  using  an  ADI  method  with  an  initial  guess  for  the  Cartesian 
velocities  and  the  pressure  field.  When  solving  these  equations,  the  contravariant  velocities 
are  first  computed  directly  from  the  guessed  Cartesian  components.  The  new  velocity  field 
will  generally  not  satisfy  the  mass  conservation  law  since  the  pressure  field  used  in  the 
solution  process  was  only  approximately  correct.  Consequently  the  pressure  field  must  be 
corrected  to  more  closely  satisfy  the  continuity  equation. 

In  the  basic  pressure  correction  procedure  of  Patankar  (1980),  the  correct  pressure 
field  p  is  obtained  from 
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p    =  p*  +p  (3.28) 

where  p*  is  an  approximate  pressure  field  and  p  is  called  the  pressure  correction. 
Corresponding  contravariant  velocity  corrections  may  be  introduced  in  a  similar  way 
following  the  procedure  described  in  Shyy  (1994)  and  Braaten  and  Shyy  (1986). 

V1   =  V\  +  V\  (3.29) 

V2    =   V*2  +  V2  (3.30) 

To  derive  the  pressure  correction  equation  we  first  subtract  the  momentum  equations  with 
the  approximate  pressure  and  velocity  fields  from  the  same  momentum  equations  using  the 
correct  flow  field  variables.  The  resulting  momentum  correction  equations  are  then  assumed 
to  be  adequately  approximated  by  the  following  truncated  form 

v\p  AlP  ~  -  x22  p,\    +  x2l  p,2  (3.31) 

v'ip  A2P  ~  xh2  p,\    -  xu  p,'2  (3-32) 

Using  these  approximate  correction  formulas,  the  Cartesian  velocity  components  may  be 
corrected  by  the  following  relationships 


v,     =    Vi    + 


(       *2,2  P'l  ~^  x2,\  P'2'  (2  23) 


Axp 


v2    =  v2  + 


(x\,2  P'\  "  x\.i  P^l)  p  34) 


2    ~    V2  A 

Subsequently,  the  contravariant  component  correction  formulas  may  be  obtained  by 
substituting  Eqs.  (3.33)  and  (3.34)  into  Eqs  (2.15)  and  (2.16).  After  dropping  terms  that  are 
not  representable  on  a  five  point  finite-difference  molecule  the  following  simplified 
correction  equations  for  the  contravariant  velocity  components  are  obtained 
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v,  =  vT  - 


'(*22)2 


+ 


i/» 


(^..2)2N 

A2P 


/M 


(3.35) 


^2  =  v5  - 


V2>1)2 


ip 


+ 


(*i.i)3 


i2P 


P.2 


(3.36) 


These  equations  are  then  substituted  into  Eq.  (3.5)  to  obtain  Eq.  (3.21). 

The  following  residuals  are  defined  to  assess  the  degree  of  convergence  of  the 
discrete  form  of  the  fluid  dynamic  conservation  laws. 


R* 


=  1 


allP 


~  A\PV\P+    X  Aln6vl«6  +  5v, 


allnb 


Qv\,c 


(3.37) 


R, 


=  1 


allP 


A2pV2p  +    2_,  A2nbv2nb  +  ^v2 


allnb 


Qv\c 


(3.38) 


R» 


=  1 


allP 


QVacC 


(3.39) 


In  the  present  procedure  the  pressure  corrections  are  used  to  update  the  contravariant 
components  using  Eqs.  (3.35)  and  (3.36).  Provided  the  pressure  correction  equation  is  solved 
to  convergence  at  each  iteration,  the  resulting  Vj  and  V?  fields  satisfy  mass  conservation 
exactly.  After  satisfying  continuity  in  terms  of  the  contravariant  components,  it  is  necessary 
to  recover  the  Cartesian  velocity  components  before  returning  to  the  momentum  equations 
with  the  updated  velocity  and  pressure  fields.  It  is  of  preeminent  importance  to  recover 
consistently  the  Cartesian  velocity  components  from  the  updated  contravariant  components 
since,  from  Eq.  (3.5),  mass  conservation  is  stated  explicitly  in  terms  of  the  contravariant 
components  in  generalized  curvilinear  coordinates.  Care  must  be  taken  to  ensure  a  consistent 
transformation  between  Cartesian  velocity  components  and  contravariant  components.  If 
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the  transformation  is  done  naively,  an  inconsistency  will  arise  in  the  discrete  form  of  the 
conservation  laws.  This  issue  is  briefly  discussed  in  the  following  paragraph. 

The  obvious  way  to  compute  the  Cartesian  components  from  the  contravariant 
components  is  to  invert  analytically  the  transformation  given  by  Eqs.  (2. 15)  and  (2.16).  This 
inversion  gives  the  following  formulas  for  the  Cartesian  components  in  terms  of  the 
contravariant  velocity  components 

X,  i  x\  2 

vi  =  viJT  +  v2~r  (3'4°) 

V2    =   y**    +   V,X~f  (3-41) 

In  order  to  estimate  the  contravariant  velocity  components,  it  is  necessary  to  interpolate  for 
the  value  of  V2  when  computing  Vj  in  Eq.  (2.15)  and  similarly  interpolate  for  the  value  of 
vi  when  computing  V2  in  Eq  (2.16).  A  corresponding  interpolation  is  needed  when 
computing  v\  and  V2  in  Eqs.  (3.40)  and  (3.41).  The  need  for  these  interpolations  causes  an 
inconsistency  in  mass  conservation  to  arise  in  the  procedure.  An  efficient  method  for 
consistently  recovering  the  Cartesian  components  from  the  updated  contravariant 
components  based  on  the  so  called  D'yakonov  iteration  procedure  is  described  in  Braaten 
and  Shyy  (1986)  and  Shyy  (1994).  In  the  present  numerical  procedure  this  method  for 
computing  the  Cartesian  velocity  components  is  adopted. 

3.3  Consistent  Implementation  of  the  Continuity  Equation 

Again,  because  of  the  need  for  interpolation  in  computing  the  contravariant  velocity 
components  appearing  in  Eqs.  (2.15-2.16)  a  possible  inconsistency  arises  in  the 
implementation  of  the  discrete  form  of  the  continuity  equation  when  the  grid  is  time 
dependent.  An  inconsistent  numerical  implementation  of  the  continuity  equation  would  lead 
to  the  generation  of  artificial  mass  sources  in  the  flow  calculation.  As  suggested  by 
Thompson  et  al.  (1985)  the  following  equation  relating  the  time  rate  of  change  of  the 
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Jacobian  to  the  Cartesian  components  of  the  grid  velocity  is  used  to  update  the  Jacobian  at 
the  implicit  time  level. 

§7  +  ;J"(~  *1*2,2  +  *2*1,2  )  +  ^f"(~  *2*1,1   +  X\X2,\    )   =  °  (3-42) 

This  identity  may  be  derived  from  Eq.  (2.14)  by  requiring  mass  to  be  conserved  for  a  constant 
density,  uniform  velocity  field  under  a  coordinate  transformation  that  is  time  dependent. 

Integrating  Eq.  (3.42)  using  a  first  order,  fully  implicit  time  integration  scheme  over 
the  same  control  volume  used  for  mass  conservation  leads  to  Eq.  (3.43)  which  is  the  discrete 
form  of  the  identity  given  above.  Adopting  this  equation  as  the  updating  formula  for  the 
Jacobian  guarantees  that  the  basic  requirement  of  geometric  conservation  (Shyy  and  Vu, 
1991,  and  Thomas  and  Lombard,  1979)  is  satisfied  in  the  discrete  form  of  the  conservation 
laws  when  the  grid  is  time  dependent. 

(J  -  7°) 

~7~  r  (  —  *i*2,2  "*"  x2x\,Ve  ~  \~  x\x2,2    '    -t2'fl,2^v*' 

+  (-  X2X\,\   +  xlx2,0n  ~  (-  *2*1,1   +  X\X2,0s  =  0  (3.43) 

3.4  The  Aeroelastic  Computational  Procedure 

The  primary  objective  of  the  present  work  is  to  determine  the  equilibrium 
configuration  and  associated  aerodynamic  characteristics  of  an  elastic  membrane  wing  as  a 
function  of  time.  Consequently,  the  present  aeroelastic  problem  consists  of  finding 
membrane  configurations  and  aerodynamic  surface  pressures  and  shear  stresses  which 
simultaneously  satisfy  Eq.  (2.1)  and  Eqs.  (2.9)  through  (2.11),  respectively.  An  iterative 
procedure  is  used  to  solve  the  coupled  boundary  value  problem  by  computing  the  elastic  and 
aerodynamic  problems  cyclically  until  a  solution  is  obtained  at  each  time  step  that  satisfies 
the  governing  equations  to  a  predetermined  convergence  criterion. 

Since  a  fully  implicit  time  integration  scheme  is  used  to  advance  the  solution  from 
one  time  level  to  the  next,  all  grid  dependent  terms  appearing  in  Eqs.  (3.3)  through  (3.5)  must 
be  reevaluated  each  time  the  membrane  profile  and  field  grid  points  are  adjusted  during  the 
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aeroelastic  iteration  procedure.  The  Cartesian  components  of  the  grid  velocity  must  also  be 
reevaluated  since  the  grid  velocity  enters  the  solution  not  only  through  the  definition  of  the 
contravariant  velocity  components  defined  by  Eqs.  (2.15)  and  (2.16),  but  also  through  the 
kinematic  boundary  condition  given  in  Eq.  (2.21).  Consequently,  the  metrics  of  the 
transformation  as  well  as  the  components  of  the  grid  velocity  vector  are  recomputed  directly 
from  the  updated  grid  coordinates  each  time  the  membrane  equilibrium  equation  is  solved. 
An  exception  to  this  strategy  is  made  when  the  Jacobian  of  the  transformation  at  the  implicit 
time  level  appearing  in  Eq.  (3.5)  is  evaluated.  This  grid  dependent  quantity  is  not  computed 
directly  from  the  updated  grid  coordinates,  but  is  computed  using  the  identity  given  by  Eq. 
(3.43). 

During  the  course  of  the  aeroelastic  iteration  procedure  it  is  necessary  to  update  the 
grid  in  the  physical  domain  so  that  it  remains  body-fitted  as  the  wing  deforms.  The 
maintenance  of  the  body-fitted  grid  is  achieved  by  updating  the  grid  points  in  the  vicinity 
of  the  membrane  during  the  course  of  iteration  according  to  the  following  formula 

xf   =  xj-"  +  g(Zj(xf  -  4'"1))  (3.44) 

where  g($2 )  is  a  general  function  that  decays  with  distance  away  from  the  membrane  surface. 
The  superscripts  appearing  in  Eq.  (3.44)  indicate  the  level  of  the  aeroelastic  iteration 
procedure.  Presently,  g(^2)  is  chosen  as  an  exponential  function  defined  as 

g(£2)  =  exp(       C[  2     j  (3.45) 

where  c  j  is  a  constant  that  depends  on  the  grid  resolution  and  the  Reynolds  number.  A  similar 
strategy  using  an  exponentially  decaying  function  for  updating  field  grid  points  while 
maintaining  a  body-fitted  curvilinear  grid  has  been  used  by  Boschitsch  and  Quackenbush 
(1993). 
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3.5  A  Potential  Flow  Model  for  Thin  Wings 

In  this  section  a  method  is  presented  for  computing  the  surface  pressure  distribution 
resulting  from  the  irrotational,  incompressible  flow  over  thin,  two-dimensional  wings  of 
arbitrary  shape.  The  reason  for  including  in  the  present  work  a  numerical  procedure  for  the 
solution  of  lifting  potential  flows  is  to  highlight  the  difference  in  aerodynamic  quantities 
predicted  by  a  viscous  flow  based  membrane  wing  model  when  compared  to  previous  work 
on  membrane  wing  aerodynamics  based  on  a  potential  description  of  the  flow  field. 

The  fundamental  assumption  concerning  the  flow  field  around  the  wing  is  that  it  is 
irrotational,  and  consequently  the  velocity  field  may  be  derived  from  a  scalar  potential. 
Implicit  in  this  assumption  is  that  the  flow  is  inviscid  and  free  of  vorticity  far  upstream.  The 
use  of  point  vortex  singularities  to  model  the  lifting  potential  flow  around  thin  wings  has  its 
historical  origin  and  justification  in  work  by  James  (1972).  A  description  of  the  method  in 
its  modern  form  may  be  found  in  Katz  and  Plotkin  (1991). 


point  vortex  panel 


Figure  3.5  Potential  flow  wing  model 

Figure  3.5  shows  a  membrane  airfoil  that  has  been  discretized  into  a  number  of  linear 
segments  or  panels.  A  typical  segment  is  composed  of  a  point  vortex  and  a  control  point 
where  the  flow  tangency  condition  is  enforced.  If  the  point  vortices  and  control  points  are 
located  at  the  local  quarter  and  three-quarter  chord  points  of  the  panel,  respectively,  the 
Kutta  condition  is  implicitly  satisfied  and  no  additional  boundary  condition  is  needed  at  the 
trailing  edge  (Katz  and  Plotkin,  1991,  James  1972).  By  modeling  the  wing  as  an  assemblage 
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of  vortex  singularity  segments  of  unknown  strength  and  enforcing  the  zero  normal  relative 
velocity  condition  at  each  control  point,  designated  as  point  i,  the  following  set  of  linear 
algebraic  equations  may  be  formed  and  solved  for  the  vortex  strength  Tat  point  j. 

A^T^  (-v„  •  n  ),  (3.46) 

In  Eq.  (3.46)  Al}  is  a  matrix  of  vortex  influence  coefficients,  defined  as  the  normal  velocity 
induced  at  control  point  i  by  a  unit  strength  vortex  at  point  j,  v^  is  the  freestream  velocity 
vector,  and  n  is  the  unit  normal  vector  at  the  control  point. 

Once  the  vortex  strengths  have  been  determined  from  Eq.  (3.46)  the  perturbation 
velocities  on  the  upper  and  lower  surfaces,  v+  and  v~  respectively,  are  given  by 


2  / 


-  n. 


(3.47) 


v      =    — 


2  / 


-  n. 


(3.48) 


where  /  is  the  length  of  the  segment  or  panel  and  nj  and  n.2  are  the  Cartesian  components  of 
the  unit  normal  vector. 

With  the  velocity  field  determined  on  the  upper  and  lower  surface  of  the  wing  the 
pressure  on  the  wing  surfaces,  p+  and  p~,  may  be  computed  directly  from  the  Bernoulli 
equation  as  follows 


±         1/2 

1      =  -kg  (  vt,  -  v 


±2 


) 


(3.49) 


The  total  fluid  velocity  on  the  upper  and  lower  wing  surfaces  is  the  sum  of  the  perturbation 
velocity  and  the  free  stream  velocity  and  is  given  by 

„±    =   Voo  +  v±  (3.50) 

It  should  be  pointed  out  that  the  leading  edge  suction  force  arising  from  the  pressure 
singularity  at  the  leading  edge  of  the  infinitely  thin  membrane  is  not  included  in  the  above 
potential  flow  model. 


CHAPTER  4 
ELEMENTARY  TEST  CASES 

4, 1  Rotated  Channel  Flow 

Before  applying  the  present  method  to  an  aeroelastic  problem,  the  performance, 
accuracy  and  convergence  properties  of  the  numerical  formulations  presented  in  Chapter  3 
are  investigated  by  applying  the  membrane  equilibrium  and  fluid  dynamic  formulations 
independently  to  certain  relevant  test  problems.  The  first  test  case  to  be  presented  is  that  of 
steady,  developing  parallel  flow  in  a  channel  that  is  rotated  45°  to  the  Cartesian  axes.  Once 
the  flow  becomes  fully  developed  the  convection  terms  vanish  in  the  Navier-Stokes 
equations  and  an  exact  solution  to  the  conservation  laws  may  be  obtained  (Schlichting  1 979). 


Figure  4.1.  Pressure  contours  (a),  and  convergence  path  (b),  for  a  45°  rotated  channel  at 
Rn=40.  The  results  shown  are  for  a  41  x  41  uniform  grid.  D'yakonov  iteration  was  used 
in  the  computation. 

Figure  4.1  (a)  shows  the  constant  pressure  contours  computed  using  the  numerical 
algorithm  described  in  the  previous  chapter  with  the  D'yakonov  iterative  procedure  used  to 
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recover  the  Cartesian  velocity  components  from  the  contravariant  components.  A  uniform 
velocity  field  is  specified  at  the  channel  inlet  and  a  0th  order  extrapolation  of  the  Cartesian 
velocity  components  is  used  at  the  channel  exit  to  enforce  a  gradient  free  boundary  condition 
in  the  streamwise  direction.  As  may  be  seen  from  the  figure  the  streamwise  pressure  gradient 
becomes  a  constant  after  some  initial  adjustment  of  the  velocity  and  pressure  fields  near  the 
inlet.  A  comparison  of  the  computed  streamwise  pressure  gradient  with  the  analytic  solution 
for  fully  developed  parallel  flow  is  shown  in  Fig.  4.3.  It  may  be  seen  from  the  figure  that  once 
the  flow  becomes  fully  developed  the  computed  pressure  gradient  is  indistinguishable  from 
the  analytic  solution.  The  residuals  of  the  discrete  form  of  the  conservation  laws,  defined 
by  Eqs.  (3.37)  through  (3.39),  are  shown  in  Fig.  4. 1  (b).  The  momentum  equations  converge 
to  a  terminal  level  of  5  x  10"4  and  the  continuity  equation  converges  to  a  terminal  level  of 
5  x  10-5.  These  levels  of  residuals  are  consistent  with  the  single  precision  floating  point 
accuracy  of  the  arithmetic  used  in  the  calculation. 
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Figure  4.2.  Pressure  contours  (a),  and  convergence  path  (b),  for  a  45°  rotated  channel  at 
Rn  =  40.  The  results  shown  are  for  a  41  X  41  uniform  grid.  D'yakonov  iteration  was 
not  used  in  the  computation. 
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Figure  4.2  shows  the  results  for  the  same  problem  but  now  the  D'yakonov  procedure 
is  not  used  and  the  Cartesian  velocity  components  are  computed  directly  from  Eqs.  (3.40) 
and  (3.41).  Interestingly,  the  pressure  contours  look  qualitatively  correct  but  the  momentum 
equation  residuals  have  converged  to  a  terminal  level  which  is  order  one.  The  reason  for  the 
dissatisfaction  of  the  momentum  equations  is  evident  in  Fig.  4.3  where  it  may  be  seen  in  that 
the  pressure  gradient  is  computed  incorrectly  when  the  D'yakonov  iterative  procedure  is  not 
used.  The  convergence  path  shown  in  Fig.  4.2  (b)  clearly  illustrates  the  inconsistency  that 
arises  in  the  transformation  due  to  the  noncollocated  velocity  components.  Fortunately,  the 
convergence  path  shown  in  Fig.  4. 1  (b)  shows  that  the  inconsistency  can  be  resolved  using 
the  D'yakonov  procedure  during  the  course  of  the  computation. 
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Figure  4.3.  Centerline  pressure  coefficient  for  the  45°  rotated  channel  at  Rn  =  40. 


4.2  Uniform  Flow  Using  a  Moving  Grid 

The  second  elementary  test  case  to  be  considered  is  the  computation  of  a  uniform 
flow  on  a  grid  that  moves  arbitrarily  in  space  as  a  function  of  time.  This  is  an  important  test 
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case  since  any  formulation  that  is  to  be  applied  to  unsteady  flow  problems  using  time 
dependent,  body-fitted  coordinates  must  admit  a  uniform  flow  field  as  a  solution  identically 
when  the  grid  coordinates  are  an  arbitrary  function  of  time. 

Figure  4.4  shows  the  grid  and  pressure  contours  after  one  time  step  for  a  uniform  flow 
field  inclined  at  45°  to  the  Cartesian  coordinate  axes.  The  initial  grid  at  time  t  is  shown  in 
Fig.  4.4  (a)  and  the  grid  at  time  t  +  At  is  shown  in  Fig.  4.4  (b).  The  boundary  conditions 
imposed  on  the  solution  at  all  computational  boundaries  are  Dirichlet  conditions  on  the 
Cartesian  velocity  components.  The  pressure  contours  computed  using  Eq.  (3.43)  to  update 
the  Jacobian  of  the  transformation  at  the  t  +  At  time  level  are  shown  in  Fig.  4.4  (c).  The 
computed  pressure  gradient  in  the  flow  is  approximately  10-7  everywhere  in  the  domain  and 
the  contours  shown  in  the  figure  result  from  roundoff  error  at  the  limit  of  single  precision 
floating  point  arithmetic.  Contrastingly,  the  pressure  contours  computed  using  Eq.  (2.20)  to 
calculate  the  Jacobian  at  the  t  +  A  t  time  level  are  shown  in  Fig.  4.4  (d).  The  pressure  gradient 
in  this  case  is  order  one,  which  clearly  demonstrates  the  inconsistency  that  arises  in  the 
computation  by  simply  taking  a  snapshot  of  the  grid  at  the  implicit  time  level  and  computing 
the  Jacobian  using  Eq.  (2.20).  Adopting  Eq.  (3.43)  as  the  updating  formula  for  the  Jacobian 
guarantees  the  basic  conservation  laws  are  satisfied  for  unsteady  flows  using  time  dependent, 
body-fitted  coordinates. 

4.3  Elastic  Membrane  Under  a  Uniform  Pressure  Load 

The  last  elementary  test  problem  to  be  investigated  is  the  deformation  of  an  elastic 
membrane  under  a  prescribed  uniform  pressure  load.  Figure  4.5  (a)  shows  the  computed 
equilibrium  shape  of  an  initially  flat  membrane  under  uniform  pressure  loading.  For  this  test 
case  the  pretension  is  chosen  to  be  zero,  consequently  the  elastic  response  of  the  membrane 
is  completely  described  by  the  nondimensional  elastic  parameter  77/ .  For  the  value  of  77/ 
shown  in  the  figure  the  numerical  algorithm  converged  to  machine  accuracy  in  less  than  100 
iterations. 
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pressure  gradient 
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Figure  4.4.  Comparison  of  Jacobian  updating  strategies  for  the  computation  of  a  45° 
uniform  flow  on  a  time  dependent  grid.  Initial  and  t=t+At  grids  are  shown  in  (a)  and 
(b),  respectively.  Computed  pressure  contours  using  Eq.  (4.43)  are  shown  in  (c),  and 

using  Eq.  (2.20)  are  shown  in  (d). 
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As  the  inextensible  limit  is  approached  and  77/  tends  to  infinity  the  geometric 
nonlinearity  intrinsic  in  the  problem  becomes  more  pronounced  and  the  elastic  boundary 
value  problem  becomes  algorithmic  ally  stiff.  Consequently,  the  use  of  underrelaxation 
becomes  essential  to  ensure  convergence  of  the  numerical  scheme.  This  is  significant  since 
most  of  the  practical  applications  for  membrane  wings,  as  well  as  the  available  experimental 
data,  are  for  flexible  but  nearly  inextensible  membranes.  In  the  aeroelastic  membrane  wing 
cases  presented  in  Chapters  5  and  6,  relaxation  factors  as  small  as  10"4  were  necessary  to 
ensure  algorithmic  stability. 
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Figure  4.5.  Computed  equilibrium  configuration,  (a),  and  midpoint  coordinate,  (b),  for 
an  initially  flat,  uniformly  loaded  elastic  membrane.  Computed  results  used  100  grid 

points  along  the  membrane  chord. 

Figure  4.5  (b)  shows  a  comparison  of  the  computed  midpoint  coordinate  of  the 
membrane  with  an  analytic  solution  for  an  elastic  membrane  given  by  Seide  (1977).  Here, 
uniform  pressure  has  been  substituted  for  the  stagnation  pressure  in  the  definition  of  111 .  The 
computed  solutions  and  the  analytic  results  are  essentially  indistinguishable. 


CHAPTER  5 
MEMBRANE  WINGS  IN  STEADY  FLOW 

5. 1  Rigid  Wing  in  a  Viscous  Fluid 

In  addition  to  the  elementary  test  cases  discussed  in  the  previous  chapter  several  other 
aspects  of  the  computational  procedure  need  to  investigated  before  the  algorithm  is  applied 
to  a  membrane  wing  problem.  Specifically,  the  effect  of  locating  the  outer  computational 
boundary  and  the  effect  of  grid  refinement  on  the  steady  state  solution  is  investigated  in  this 
section. 

5.1.1  Effect  of  Outer  Boundary  Location 

Since  the  physical  problem  under  consideration  is  the  flow  over  a  wing  in  an 
unbounded  domain  the  placement  and  boundary  conditions  imposed  on  the  outer  boundary 
of  the  computational  domain  require  investigation.  Consequently,  a  sequence  of 
computations  was  performed  to  assess  the  sensitivity  of  the  computed  results  to  the  location 
and  the  boundary  conditions  imposed  at  the  outer  flow  boundary. 

Figure  5.1  (a)  shows  a  typical  grid  used  in  the  membrane  wing  computations.  The 
wing  is  situated  in  the  center  of  a  square  computational  domain  typically  10c  x  10c  in  size. 
The  freestream  velocity  is  imposed  on  all  outer  flow  boundaries  except  the  downstream 
boundary  where  a  zero  gradient  condition  is  imposed  on  the  velocity  components.  An 
enlarged  view  of  the  grid  in  the  vicinity  of  the  wing  is  shown  in  Fig  5. 1  (b).  The  refinement 
of  the  grid  near  the  membrane  was  based  primarily  on  experience  and  the  degree  to  which 
flow  separation  was  expected  in  each  situation.  Attention  to  grid  refinement  is  particularly 
important  near  the  trailing  edge  and  in  the  wake  where  large  regions  of  recirculating  flow 
often  occurred. 
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Figure  5.1.  Typical  grid  used  for  membrane  wing  computations.  The  complete 
computational  domain  and  boundary  conditions  are  shown  in  (a),  and  an  enlarged  view 
of  the  grid  near  the  membrane  is  shown  in  (b). 

A  series  of  computations  was  performed  using  different  conditions  at  the  outer 
domain  boundary  prior  to  selecting  the  boundary  condition  set  shown  in  Fig.  5.1  (a).  These 
included  computations  using  Dirichlet  conditions  at  all  outer  domain  boundaries  as  well  as 
computations  where  a  zero  gradient  condition  was  imposed  on  all  boundaries  except  the 
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upstream  boundary  where  freestream  conditions  were  imposed.  In  the  latter  case  with  exit 

conditions  specified  at  the  north,  south,  and  east  domain  boundaries,  the  algorithm  was 

unable  to  admit  a  uniform  flow  field  as  a  solution  to  the  fluid  dynamic  conservation  laws. 

Consequently,  the  exit  boundary  condition  was  retained  only  at  the  downstream  domain 

boundary  as  shown  in  Fig.  5.1  (a). 

Table  5.1.  Effect  of  outer  boundary  location  on  computed  aerodynamic  properties  of  a 
rigid,  2%  camber,  zero  thickness,  circular  arc  airfoil  at  a  =  5°,  Rn  =  4xl03. 


grid 

domain  size 

CL 

CD 

Cm 

185X77 

5cx5c 

.637 

.0761 

.0442 

207x91 

10c X 10c 

.626 

.0752 

.0425 

221x101 

15cXl5c 

.623 

.0750 

.0422 

247x121 

30c X 30c 

.621 

.0749 

.0420 

The  effect  of  moving  the  outer  boundary  away  from  a  rigid,  zero  thickness  flat  plate 
at  4°  angle  of  attack  and  a  Reynolds  number  of  4x  1 03  may  be  seen  in  Table  5.1.  The  boundary 
was  moved  outward  by  adding  additional  grid  points  to  the  basic  5c  X  5c  grid  so  that  the  outer 
boundary  location  is  isolated  from  other  grid  dependency  issues.  For  all  grids  100  grid  points 
were  used  along  the  wing  chord.  As  may  be  seen  from  the  table  there  is  very  little  change 
-  less  than  1%  -  in  the  computed  lift,  drag  and  aerodynamic  moment  beyond  a  domain  size 
of  10c  X  10c.  Consequently,  a  computational  domain  size  of  10c  X  10c  with  exit  conditions 
imposed  at  the  downstream  boundary  is  adopted  for  all  further  computations. 

5.1.2  Effect  of  Grid  Refinement 


In  order  to  asses  the  effect  of  grid  resolution  on  the  computed  flow  field  near  a 
membrane  wing  two  test  cases  were  chosen  for  investigation.  The  first  test  case  was  a  rigid, 
zero  thickness  flat  plate  at  4°  angle  of  attack  and  a  Reynolds  number  of  4xl03.  Figure  5.2 
shows  the  computed  streamlines  for  the  flat  plate  using  100,  200,  and  300  grid  points  along 
the  wing  chord.  As  may  be  seen  in  the  figure  a  small  separation  bubble  emerges  at  the  leading 
edge  and  grows  as  the  grid  is  refined.  The  computed  lift,  drag,  and  moment  taken  about  the 
wing  quarter  chord  are  presented  in  Table  5.2  for  the  grid  refinement  sequence  of  Fig  5.2. 
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Table  5.2.  Effect  of  grid  refinement  on  computed  aerodynamic  properties  of  a  rigid  flat 

plate  at  a  =  4°,  7?/z  =  4xl63 

grid 

points  on  c 

CL 

CD 

Cm 

161x81 

100 

.401 

.0673 

.0012 

281x101 

200 

.404 

.0604 

.0009 

401  X 121 

300 

.413 

.0572 

.0022 

The  table  shows  that  drag,  as  expected,  is  the  quantity  most  strongly  affected  by  grid 
resolution.  The  computed  lift  varies  by  less  than  3%  over  the  range  of  grid  resolution 
investigated.  The  computed  aerodynamic  moment  is  effectively  zero  for  all  grids.  This  is 
consistent  with  the  analytic  result  from  thin  wing  theory  (Katz  and  Plotkin  1991). 


Figure  5.2.  Effect  of  grid  refinement  on  streamfunction  contours  for  a  rigid  flat  plate  at 
a=4°, /?/z=4xl03. 
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Figure  5.3  shows  the  computed  aerodynamic  surface  pressures  at  Rn  =  104  for  a  rigid 
flat  plate  at  3°  angle  of  attack.  The  computed  surface  pressures  used  a  grid  with  100  points 
along  the  wing  chord.  Also  shown  in  the  figure  for  comparison  is  the  potential  flow  solution 
obtained  from  conformal  mapping  for  the  flat  plate.  The  analytical  solutions  for  the  potential 
flow  is  taken  from  Katz  and  Plotkin  ( 199 1 ).  The  computed  surface  pressures  for  the  flat  plate 
are  in  good  agreement  with  the  potential  flow  solution  since  the  flow  is  largely  attached  and 
the  boundary  layer  is  fairly  thin  at  this  Reynolds  number. 
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Figure  5.3.  Comparison  of  computed  surface  pressures  at  Rn=l04  and  a=3°  with 
potential  flow  solution  for  a  rigid,  zero  thickness  flat  plate. 


Figure  5.4  shows  the  computed  lift  curves  at  Reynolds  numbers  of  104  and  4  x  103 
for  a  rigid  flat  plate  using  100  grid  points  along  the  wing  chord.  Also  shown  in  the  figure 
for  comparison  is  the  potential  flow  lift  curve  taken  from  Katz  and  Plotkin  (1991).  Again, 
as  in  Fig.  5.3,  the  viscous  flow  solution  shows  close  agreement  with  the  potential  flow 
solution  at  these  Reynolds  numbers  and  angles  of  attack. 
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Figure  5.4.  Comparison  of  computed  lift  coefficient  with  potential  flow  solution  for  a 
rigid,  zero  thickness  flat  plate. 


The  second  test  case  for  investigating  the  effect  of  grid  refinement  on  the  computed 
flow  field  is  a  rigid,  2%  camber,  circular  arc  airfoil  at  5°  angle  of  attack  and  a  Reynolds 
number  of  4x  1 03.  Figure  5.5  shows  the  computed  streamlines  for  the  circular  arc  airfoil  using 
100,  200,  300  and  400  grid  points  along  the  wing  chord.  As  was  seen  with  the  flat  plate  grid 
refinement  sequence  a  leading  edge  separation  bubble  appears  as  the  grid  is  refined. 
However,  the  most  significant  flow  feature  to  emerge  as  the  grid  is  refined  is  the  recirculation 
region  near  the  trailing  edge  and  the  attendant  departure  of  the  streamlines  from  the  upper 
surface  of  the  airfoil.  The  wake  behind  the  airfoil  may  also  be  seen  to  thicken  significantly 
with  the  appearance  of  the  trailing  edge  separation  region.  Based  on  the  grid  refinement 
sequence  shown  in  Figure  5.5,  it  may  be  concluded  that  for  circular  arc  configurations  - 
associated  with  either  flexible  or  rigid  wings  -  potential  flow  will  not  be  capable  of 
accurately  describing  the  flow  field. 
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Figure  5.5.  Effect  of  grid  refinement  on  streamfunction  contours  for  a  rigid,  2%  circular 
arc  airfoil  at  a=  5°,  Rn=4x\03. 


The  computed  lift,  drag,  and  aerodynamic  moment  for  the  circular  arc  airfoil  are 
presented  in  Table  5.3  for  the  grid  refinement  sequence  shown  in  Fig  5.5.  In  contrast  to  the 
flat  plate  results  the  computed  values  for  the  lift,  drag,  and  aerodynamic  moment  are  quite 
sensitive  to  grid  refinement.  Even  at  the  highest  grid  resolution  there  is  a  5%  -  10%  change 
in  computed  values  when  compared  to  the  previous  grid.  However,  it  is  expected  that  a  grid 
independent  solution  will  be  difficult  to  obtain  for  this  problem  since  a  stringent  requirement 
is  placed  on  grid  resolution  at  this  Reynolds  number. 
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Table  5.3.  Effect  of  grid  refinement  on  computed  aeroelastic  properties  for  a  rigid,  2% 

circular  arc  airfoil  at  a  =  5°,  Rn  =  4xl03 


grid 

points  on  c 

CL 

CD 

Cm 

207x91 

100 

.626 

.0752 

.0425 

311x93 

200 

.595 

.0654 

.0394 

411x95 

300 

.567 

.0603 

.0378 

531x105 

400 

.531 

.0558 

.0352 

Based  on  the  results  for  a  circular  arc  airfoil  it  may  be  expected  that  computations 
with  a  flexible  membrane  will  exhibit  a  similar  sensitivity  to  grid  refinement  since  the 
equilibrium  configuration  for  a  pressure  loaded,  end  restrained  membrane  is  quite  similar 
to  a  circular  arc  and  will  likely  exhibit  the  same  type  of  trailing  edge  separation. 
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Figure  5.6.  Comparison  of  computed  surface  pressures  at  Rn=l04  and  a=3°  with 
potential  flow  solution  for  a  rigid,  zero  thickness,  2%  camber,  circular  arc  airfoil. 

Figure  5.6  shows  the  computed  aerodynamic  surface  pressures  at  Rn  =  104  for  a  rigid, 
2%  camber  circular  airfoil  at  3°  angle  of  attack.  The  computed  surface  pressures  used  a  grid 
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with  100  points  along  the  wing  chord.  Again,  the  potential  flow  solution  obtained  from 
conformal  mapping  for  a  circular  arc  airfoil  taken  from  Katz  and  Plotkin  (1991)  is  shown 
in  the  figure  for  comparison.  The  computed  surface  pressures  for  the  arc  airfoil  are  in 
reasonable  agreement  with  the  potential  flow  solutions.  Figure  5.7  shows  the  computed  lift 
curves  using  100  grid  points  along  the  wing  chord  at  Reynolds  numbers  of  104  and  4xl03 
for  the  rigid,  circular  arc  airfoil.  Inspection  of  Figs.  5.4  and  5.7  indicate  there  is  a  general 
loss  of  circulation,  compared  to  potential  flow  theory,  around  the  airfoil  due  to  viscous 
effects  near  the  trailing  edge.  The  influence  of  viscosity  is  much  more  pronounced  for  the 
circular  arc  than  for  the  flat  plate  at  this  angle  of  attack  and  Reynolds  number. 
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Figure  5.7.  Comparison  of  computed  lift  coefficient  with  potential  flow  solution  for  a 
rigid,  zero  thickness,  2%  camber  circular  arc  airfoil. 


5.2  Flexible  Membrane  Wing  in  a  Viscous  Fluid 

Membrane  wings  may  be  broadly  classified  by  considering  the  following  three 
limiting  cases  of  the  parameters  77/,  772,  and  e.  The  first  limiting  case  may  be  referred  to 
as  the  elastic  wing  case  with 
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77/  _  finite 

n2  =0 

e  =0 
In  this  case  the  wing  is  initially  flat  and  taut  but  with  no  pretension.  Consequently,  the 
membrane  behavior  is  determined  exclusively  by  77/ . 

The  second  limiting  case  for  membrane  wings  may  be  referred  to  as  the  constant 
tension  wing  case  with 

77/  =0 

TTj  -  finite 

£=0 

As  with  the  first  case  the  wing  is  initially  flat  and  taut  but  with  now  with  sufficient  pretension 
so  that  the  material  stiffness  does  not  participate  in  determining  the  equilibrium  membrane 
configuration. 

The  third  limiting  case  may  be  referred  to  as  the  inextensible  wing  case  with 

77;  _  infinite 

772  =0 

e  _  finite 
In  this  limiting  case  the  wing  is  not  initially  flat  and  taut  but  has  an  excess  length  of 
material  defined  by  the  parameter  e.  Of  course,  the  Reynolds  number  and  the  angle  of  attack 
are  additional  parameters  that  characterize  a  steady  flow  about  an  airfoil,  but  they  are  to  be 
distinguished  from  the  three  parameters  listed  above  that  characterize  flexible  membrane 
wing  behavior. 

5.2.1  Elastic  Membrane  Case 

The  computed  streamlines  and  constant  pressure  contours  for  an  elastic  membrane 
wing  typical  of  the  first  limiting  case  with  Tli=10, 112=0,  £=0,  are  shown  in  Fig.  5.8.  The 
solution  shown  in  the  figure  was  computed  using  a  grid  with  100  points  along  the  wing  chord. 
The  recirculation  region  near  the  trailing  edge  is  reminiscent  of  the  flow  over  a  rigid  circular 
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arc  airfoil  as  discussed  above.  Consequently,  it  may  be  anticipated  that  further  refinement 
of  the  grid  would  lead  to  the  emergence  of  additional  flow  features,  such  as  a  leading  edge 
separation  bubble,  as  was  demonstrated  with  the  rigid  arc  airfoil. 


Figure  5.8.  Streamlines  and  constant  pressure  contours  for  an  initially  flat  elastic 
membrane  wing,  Rn=4x\02,  a=6°  777=10, 772=0,  e=0. 


Figure  5.9  illustrates  the  convergence  path  of  the  four  coupled  governing  equations 
for  the  the  elastic  wing  case  of  Fig.  5.8.  It  may  be  seen  that  the  residuals,  as  defined  by  Eq. 
(3.2)  and  Eqs.  (3.37)  through  (3.39)  are  reduced  to  single  precision  machine  accuracy  after 
a  few  thousand  iterations.  Also,  the  terminal  level  of  the  residuals  shown  in  the  figure  is 
consistent  with  the  single  precision  floating  point  accuracy  of  the  arithmetic  used  in  the 
calculation.  It  may  be  noted  that  the  stability  of  the  aeroelastic  algorithm  is  largely  dictated 
by  the  choice  of  relaxation  parameter  used  in  the  solution  of  the  membrane  equilibrium 
equation.  For  the  solution  shown  in  the  Fig.  5.8  and  Fig.  5.9,  a  relaxation  value  of  10-2  was 
used  during  the  solution  of  Eq.  (3.1). 
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Figure  5.9.  Convergence  path  of  the  momentum,  continuity,  and  equilibrium  equations 
for  an  initially  flat  elastic  membrane  wing,  /?/i=4xl03,  a=6°,  77/  =10,  772=0,  e=0. 


Figure  5. 10  compares  the  computed  membrane  surface  pressures  for  the  case  of  Fig. 
5.8  with  the  surface  pressures  predicted  by  a  potential  flow  calculation  for  the  same 
membrane  configuration.  Again,  as  with  Fig.  5.7,  the  discrepancy  between  the  potential 
solution  and  viscous  flow  solutions  appears  to  be  attributable  to  a  general  loss  of  circulation 
about  the  wing  due  to  viscous  effects  primarily  near  the  trailing  edge.  However,  for  this  set 
of  aeroelastic  parameters,  potential  flow  does  reasonably  well  in  predicting  the  membrane 
surface  pressures.  For  other  choices  of  the  dimensionless  parameters,  potential  flow  is  not 
as  successful  in  accurately  describing  the  fluid  dynamics  of  the  problem.  In  the  following 
section  describing  the  aeroelastic  behavior  of  an  inextensible  membrane  wing  of  moderate 
excess  length,  the  limitations  of  a  potential  based  description  of  the  flow  field  will  be 
demonstrated. 
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Figure  5.10.  Equilibrium  surface  pressures  for  an  initially  flat  elastic  membrane  wing, 
Rn=4xl03,  a=6°,  77/=10, 772=0,  6=0.  Also  shown  is  the  potential  flow  surface  pressure 
computed  for  the  same  membrane  configuration. 


Figure  5. 1 1  (a)  shows  the  equilibrium  lift  coefficient  for  the  initially  flat  elastic  wing 
of  Fig.  5.8.  The  figure  shows  the  dependency  on  Reynolds  number  to  be  more  pronounced 
for  the  elastic  wing  when  compared  to  the  rigid  circular  arc  airfoil  of  Figure  5.7.  This  may 
be  explained  by  observing  that  the  maximum  wing  camber,  shown  in  Fig.  5.11  (b),  follows 
the  same  trend  as  the  lift  coefficient.  It  may  be  concluded  that  this  enhanced  Reynolds 
number  dependency  is  largely  a  result  of  the  coupling  between  the  membrane  deformation 
and  the  flow  field.  Also  shown  in  the  figure  is  the  equilibrium  lift  coefficient  predicted  by 
a  potential  based  membrane  wing  model  for  the  same  set  of  aeroelastic  parameters.  It  may 
be  seen  that  the  potential  based  solution  considerably  overpredicts  the  lift  on  the  wing  even 
at  small  camber  ratios  and  angles  of  attack  when  compared  to  the  viscous  flow  based  solution 
at  these  Reynolds  numbers.  As  stated  previously  this  discrepancy  can  be  attributed  to  a 
general  loss  of  circulation  about  the  wing  due  to  viscous  effects  near  the  trailing  edge. 
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Figure  5.11.  Equilibrium  lift  coefficient,  (a),  and  X2  coordinate  at  midchord,  (b),  for  an 
initially  flat  elastic  membrane  wing,  77; =10, 772=0,  e=0. 


5.2.2  Constant  Tension  Membrane  Case 

Figure  5. 12  shows  the  equilibrium  lift  coefficient  and  maximum  wing  camber  for  the 
second  limiting  case  of  a  constant  tension  membrane  with  77;  =0, 772=2,  and  e=0.  Again,  1 00 
grid  points  were  used  along  the  membrane  chord.  This  figure  may  be  contrasted  with  Fig. 
5.11.  The  contrast  between  constant  tension  and  elastic  wing  characteristics  may  be 
reconciled  by  recalling  that  the  response  of  an  initially  flat,  constant  tension  membrane  to 
a  specified  uniform  pressure  load  is  qualitatively  different  from  the  response  of  a  membrane 
with  finite  material  stiffness.  For  small  deflections  the  constant  tension  membrane  will 
exhibit  linear  elastic  behavior  whereas  an  elastic  membrane  with  finite  material  stiffness  will 
exhibit  a  1/3  power  law  response  to  the  applied  pressure  load  (Seide  1977)  as  shown  in  Fig. 
4.5.  An  inspection  of  Figs.  5.12  and  5.11  illustrates  the  different  behavior  of  the  two 
membrane  wings  near  zero  degrees  incidence.  The  initially  flat,  unpretensioned  elastic 
membrane  wing  of  Fig.  5.11  has  very  little  initial  resistance  to  deformation,  and 
consequently,  a  substantial  amount  of  camber  exists  in  the  equilibrium  configuration  near 
zero  degrees  angle  of  attack.  Conversely,  the  pretensioned  membrane  wing  of  Fig.  5.12  has 
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adequate  stiffness  in  its  initially  flat  unstrained  configuration  to  postpone  the  development 
of  camber  until  higher  angles  of  attack  are  reached.  For  the  constant  tension  wing  of  Fig.  5.12 
the  lift  curve  closely  mimics  the  the  development  of  camber  with  angle  attack  as  was  the  case 
with  the  elastic  wing  of  Fig.  5.11.  Streamlines  and  pressure  contours  are  not  presented  for 
this  limiting  case  since  they  are  similar  to  the  first  limiting  case  shown  in  Fig.  5.12. 


0.08 


Figure  5.12.  Equilibrium  lift  coefficient,  (a),  and  xj  coordinate  at  mid  chord,  (b),  for  an 
initially  flat,  constant  tension  membrane  wing,  77;  =0, 112=2,  e=0. 

5.2.3  Inextensible  Membrane  Case 

Finally,  the  limiting  case  of  a  flexible  but  inextensible  membrane  wing  with  excess 
length  e=. 017, 77/ =46  and 772=0  is  illustrated  in  Fig.  5.13.  The  excess  length  and  aeroelastic 
parameters  were  chosen  to  reproduce  a  portion  of  the  experimental  data  reported  by  Newman 
and  Low  (1984).  Experience  has  shown  that  for  this  excess  length  any  further  increase  in  77/ 
will  not  substantially  effect  the  aeroelastic  solution.  Consequently,  the  results  shown  in  Fig. 
5.13  are  for  a  membrane  wing  that  may  be  considered  essentially  inextensible.  Due  to  the 
assumption  of  steady  laminar  flow  the  numerical  computations  were  limited  to  Reynolds 
numbers  below  4xl03  whereas  the  Reynolds  number  reported  by  Newman  &  Low  (1984) 
was  1.2xl05  which  is  in  the  turbulent  flow  regime. 
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300  points  on  c 


400  points  on  c 


Figure  5.13.  Effect  of  grid  refinement  on  streamfunction  contours  for  an  inextensible 
membrane  wing  at  a=0°,  77/=46,  772=0,  e=.017,  Rn=4\l03. 


The  effect  of  grid  refinement  on  the  computed  streamlines  for  the  inextensible 
membrane  wing  at  Rn=4x\03  and  0°  angle  of  attack  is  shown  in  Fig.  5.13.  As  may  be  seen 
in  the  figure  a  large  recirculation  region  appears  on  the  lower  surface  of  the  membrane  as 
the  grid  is  refined.  The  integrated  aerodynamic  properties  and  the  nominal  membrane 
tension  are  given  in  Table  5.4  for  the  grid  refinement  sequence  of  Fig.  5. 13.  Even  for  the 
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highest  grid  resolution  shown  the  tabulated  values  of  lift  and  tension  are  not  yet  grid 
independent.  Of  course,  the  membrane  lift  and  tension  will  follow  the  same  trend  as  the  grid 
is  refined  as  required  by  equilibrium.  The  highest  resolution  grid  shown  in  the  figure 
required  24  hours  of  CPU  time  on  DEC  3000  AXP  workstation. 


Table  5.4.  Effect  of  grid  refinement  on  computed  aeroelastic  properties  for  an 
inextensible  membrane  wing  at  a=0°,  77^=46,  772=0,  e=.017,  /?n=4xl03 


grid 

points  on  c 

CL 

CD 

Cm 

CT 

max  X2 

261x101 

100 

.378 

.0716 

.152 

.893 

.0836 

361x121 

200 

.342 

.0689 

.144 

.823 

.0835 

461x141 

300 

.317 

.0685 

.139 

.773 

.0835 

561x141 

400 

.300 

.0702 

.136 

.741 

.0834 

Figure  5.14  contrasts  the  viscous  flow  based  aeroelastic  solution  using  400  points 
along  the  wing  chord  with  a  potential  flow  based  aeroelastic  solution  for  the  inextensible 
membrane  wing  of  Fig.  5.13.  Although  the  equilibrium  membrane  profiles  are  similar,  as 
shown  in  Fig.  5.14  (a),  the  surface  pressures,  shown  in  Fig.  5.14  (b),  are  dramatically 
different.  In  particular  near  the  leading  edge  the  pressure  jump  across  the  membrane,  Ap, 
predicted  by  the  two  flow  models  are  of  different  sign.  The  negative  Ap  computed  by  the 
viscous  flow  model  implies  an  inflection  point  in  the  equilibrium  membrane  profile  which, 
upon  close  inspection,  may  be  seen  in  the  figure.  In  contrast  the  potential  based  model 
predicts  a  flow  pattern  and  a  membrane  profile  that  is  completely  symmetric  about  the 
midchord.  It  may  be  concluded  that  for  this  set  of  dimensionless  parameters,  a  potential  flow 
based  membrane  wing  theory  essentially  fails. 

Figure  5.15  (a)  compares  the  computed,  laminar  flow  equilibrium  lift  coefficient  at 
Reynolds  numbers  of  2xl03  and  4xl03  with  the  experimental  data  of  Newman  and  Low 
(1984)  who  reported  a  Reynolds  number  of  1.2xl05  for  their  experiments.  The  solution 
shown  in  the  figure  was  computed  using  a  grid  with  100  points  along  the  wing  chord.  The 
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disagreement  between  the  computational  and  experimental  results  is  largely  due  to  the 
assumption  of  laminar  flow  in  the  present  computation. 
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Figure  5.14.  Comparison  of  potential  and  viscous  flow  based  membrane  configurations, 
(a),  and  surface  pressures,  (b),  for  an  inextensible  membrane  wing  with  excess  length. 
The  dimensionless  parameters  are  77/ =46, 772=0,  e=.017,  a=0°. 
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Figure  5.15.  Comparison  of  lift  coefficient  with  experimental  data,  (a),  and  separation 
point,  (b),  for  an  inextensible  membrane  wing  with  excess  length,  77/=46,  772=0,  E=.017 
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The  decrease  in  lift  with  increasing  Reynolds  number  shown  in  Fig.  5.15  (a)  is 
contrary  to  the  trend  observed  in  the  previous  computations  on  membrane  wings.  This  trend 
may  be  explained  by  observing  that  the  separation  point,  xjsP  (defined  as  the  point  of 
vanishing  shear  stress),  moves  forward  along  the  membrane  as  the  Reynolds  number  is 
increased  as  shown  in  Fig.  5.15  (b).  The  presence  of  turbulent  mixing  in  the  boundary  layer 
would  undoubtedly  postpone  flow  separation  and  lead  to  higher  lift  coefficients  consistent 
with  the  experimental  data.  Also  shown  in  Fig.  5. 15  (a)  is  the  computed  equilibrium  lift  using 
a  potential  based  membrane  wing  theory.  The  potential  flow  based  aeroelastic  solution  is 
seen  to  considerably  overestimate  the  lift  on  the  membrane  wing  when  compared  to  the 
experimental  results.  Conversely,  the  viscous,  laminar  flow  based  model  considerably 
underestimates  the  lift  reported  by  Newman  and  Low  (1984). 


CHAPTER  6 
MEMBRANE  WINGS  IN  UNSTEADY  FLOW 

In  this  chapter  the  role  of  viscosity  in  the  unsteady  flow  over  a  membrane  wing  is 
investigated.  In  many  situations  viscosity  may  play  a  pronounced  role  in  the  unsteady 
aerodynamics  of  membrane  wings  since  the  membrane  will  adjust  to  changes  in  the 
freestream  velocity  as  required  by  equilibrium.  This  adjustment  in  the  membrane  shape  may 
lead  to  the  periodic  appearance  or  disappearance  of  regions  of  separated  flow.  In  addition 
to  investigating  the  role  of  viscosity  in  the  aeroelastic  boundary  value  problem,  the 
simulation  of  a  membrane  wing  subjected  to  a  harmonically  varying  freestream  velocity  also 
serves  to  demonstrate  the  differences  between  the  three  limiting  aeroelastic  cases  discussed 
in  the  previous  chapter.  On  a  more  practical  note,  the  simulation  of  a  membrane  wing 
subjected  to  a  harmonically  varying  freestream  is  a  rudimentary  yet  meaningful  model 
problem  for  investigating  the  behavior  of  a  marine  sail  in  a  wind  gust.  For  each  of  the  three 
limiting  cases  presented  in  the  following  sections  the  freestream  velocity  varies  by  20% 
about  the  mean  value. 

The  time  step  required  to  resolve  adequately  the  major  features  of  the  unsteady  flow 
about  a  membrane  wing  may  be  determined  by  comparing  the  aeroelastic  solutions 
computed  using  several  different  time  steps.  The  constant  tension  wing  was  chosen  for  this 
study  since  experience  has  shown  this  case  to  be  quite  responsive  to  a  harmonically  varying 
freestream  velocity. 

Figure  6. 1  shows  the  time  series  of  the  aeroelastic  solution  using  dimensionless  time 
steps  of  .075,  .15,  and  .30  with  a  281x101  spatial  grid.  The  freestream  velocity, 
aerodynamic  lift  coefficient,  and  membrane  X2  coordinate  at  the  midchord  point  for  each  of 
the  three  time  steps  are  shown  in  the  figure.  The  dimensionless  parameters  that  govern  the 
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solution  are  given  in  the  figure  caption.  It  may  be  seen  from  Fig.  6. 1  that  even  for  the  smallest 
time  step  the  aeroelastic  solution  is  not  yet  time  step  independent.  It  is  interesting  to  note  that 
the  solution  is  essentially  independent  of  the  time  step  during  the  portion  of  the  cycle  when 
the  freestream  is  accelerating,  and  conversely,  strongly  dependent  on  the  time  step  when  the 
freestream  is  decelerating.  The  reason  for  this  strong  timestep  dependency  during  periods 
of  freestream  deceleration  will  be  discussed  in  the  following  section.  However,  on  a  more 
practical  note,  it  is  known  from  the  previous  chapter  that  the  steady  state  solution  is  also  not 
completely  grid  independent  at  this  level  of  spatial  resolution.  Therefore,  the  spatial  and 
temporal  accuracy  afforded  by  a  281  X  101  grid  and  a  nondimensional  time  step  of  .075  is 
considered  adequate  for  the  present  investigation  and  is  adopted  for  all  unsteady 
computations.  This  choice  is  simply  an  accommodation  to  limitations  in  the  available 
computing  resources. 


i  CL,  At=.075 
-a  x^IOISc/Z,  At=.075 
-■CL,AI=15 
-»  Xj*10@c/2,  At=.  1 5 

•  q,  At=.30 
-*x*10@c/2,At=.30 
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Figure  6.1.  Variation  in  the  computed  aerodynamic  lift  and  midchord  X2  coordinate 
using  three  different  timesteps  for  a  constant  tension  membrane  wing  subjected  to 
harmonic  freestream  oscillation.  The  nondimensional  parameters  are  e=0,  /?«=4xl03, 
5f=1.5, 77/  =0,  n2=2,  and  a=4°. 
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6.1  Constant  Tension  Membrane  Case 

The  first  case  to  be  investigated  is  that  of  an  initially  flat,  constant  tension  wing. 
Figure  6.2  (a)  shows  the  time  series  of  the  free  stream  velocity,  the  aerodynamic  lift,  drag 
and  moment,  and  the  X2  coordinate  at  the  midchord  point  for  the  membrane  wing  at  4°  angle 
of  attack.  The  initial  condition  for  this  simulation  is  taken  to  be  the  steady  state  solution  to 
the  aeroelastic  problem  with  the  same  dimensionless  parameters.  Also,  the  time  shown  in 
Fig.  6.2  (a)  is  the  dimensionless  time  defined  by  Eq.  (2.26).  Figure  6.2  (b)  and  6.2  (c)  show 
the  membrane  surface  pressures  and  membrane  profile,  respectively,  at  several  instants 
during  the  second  complete  cycle  of  the  freestream  velocity. 

Several  observations  may  be  made  concerning  the  solution  shown  in  Fig.  6.2  (a). 
First,  it  may  be  seen  that  the  peak  in  the  membrane  deflection  lags  the  peak  in  the  freestream 
velocity  by  approximately  60°.  This  phase  lag,  as  well  as  the  large  amplitude  of  the  motion 
of  the  membrane  -  the  membrane  camber  varies  from  negative  1%  to  positive  8%  -  is 
characteristic  of  a  system  that  is  being  driven  near,  but  below  the  system  natural  frequency. 
This  observation  is  supported  by  the  value  of  the  frequency  ratio  Q2  (Eq.  2.49),  which  has 
a  value  that  is  approximately  1 .0  for  this  case.  Secondly,  the  membrane  profiles  shown  in 
Fig.  6.2  (c)  demonstrate  that  the  algorithm  is  capable  of  successfully  simulating  the 
dynamics  of  membrane  wings  when  substantial  changes  in  the  membrane  profile  occur. 

The  variation  in  the  membrane  profile  may  also  be  seen  in  Fig.  6.3  where  the 
streamfunction  contours  are  shown  at  equally  spaced  intervals  during  the  second  complete 
cycle  of  the  freestream  velocity.  In  this  figure  it  may  be  seen  that  as  the  freestream  velocity 
decelerates  the  flow  separates  along  the  upper  surface  of  the  membrane  and  two  regions  of 
recirculating  flow  develop  near  the  trailing  edge.  The  periodic  appearance  and  collapse  of 
these  recirculation  zones,  as  seen  in  Fig.  6.3,  suggests  the  role  of  viscosity  is  enhanced  in  the 
unsteady  flow  scenario  since  the  acceleration  and  deceleration  of  the  freestream  velocity 
strongly  influences  the  separation  and  reattachment  of  the  flow.  The  periodic  appearance  and 
collapse  of  these  flow  features,  along  with  an  attendant  adjustment  in  the  membrane 
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configuration,  results  in  an  aeroelastic  response  which  may  not  be  characterized  as  simple 
harmonic  response  at  the  freestream  forcing  frequency.  A  leading  edge  separation  bubble 
may  also  be  seen  to  appear  near  the  end  of  the  freestream  cycle  and  then  collapse  as  the 
freestream  velocity  accelerates.  The  associated  pressure  contours  are  shown  in  Fig.  6.4. 

6.2  Elastic  Membrane  Case 

The  second  unsteady  case  to  be  investigated  is  that  of  an  initially  flat  elastic  wing. 
In  this  case  the  membrane  develops  tension  as  a  result  of  material  strain  rather  than 
pretension.  The  controlling  aeroelastic  parameter  U )  was  chosen  so  that  the  steady  state 
aeroelastic  solution  is  the  same  as  the  constant  tension  case  previously  discussed.  Here,  as 
with  the  previous  case,  the  steady  state  solution  was  used  to  initialize  the  solution  at  time  t =0. 

The  time  series  of  the  aerodynamic  and  membrane  properties  are  shown  in  Fig.  6.5 
(a).  Again,  several  observations  may  be  made  concerning  this  time  series.  First,  the  lift  and 
membrane  deflection  show  very  little  deviation  from  the  steady  state  value  when  compared 
to  the  deviations  observed  in  the  constant  tension  case.  The  difference  between  the  two  cases 
may  be  explained  by  recalling  that  the  deflection  of  an  initially  flat,  elastic  membrane 
without  pretension  is  proportional  to  the  cube  root  of  the  applied  pressure.  Consequently,  an 
elastic  wing  will  become  progressively  stiffer  as  the  freestream  velocity  is  increased.  This 
effect  is  usually  referred  to  as  geometric  or  stress  stiffening  (Leonard  1990).  In  this  case  the 
effect  of  geometric  stiffening  serves  to  restrict  the  development  of  wing  camber,  and 
consequently,  the  appearance  of  regions  of  separated  flow  is  less  pronounced  than  in  the 
constant  tension  case  of  the  previous  section.  The  membrane  deflection,  as  well  as  the 
aerodynamic  properties,  may  also  be  seen  to  be  essentially  simple  harmonic  response  at  the 
freestream  forcing  frequency. 

An  interesting  feature  of  the  elastic  wing  case  is  the  apparent  lack  of  a  phase  shift 
between  the  freestream  velocity  and  the  membrane  deflection.  This  situation  is  consistent 
with  a  system  that  is  being  driven  well  below  its  natural  frequency.  This  observation  is 
supported  by  the  frequency  ratio  Q\  (Eq.  2.48),  which  for  this  problem  is  approximately  0. 1 . 
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Figure  6.2.  Time  series  of  system  response  for  a  constant  tension  membrane  wing 
subjected  to  harmonic  freestream  oscillation  is  shown  in  (a).  Surface  pressures  and 
membrane  profiles  at  several  times  during  one  complete  oscillation  of  the  freestream 
velocity  are  shown  in  (b)  and  (c),  respectively.  The  nondimensional  parameters  are 
e=0,  toz=4xl03,  Sr=1.5,  Ui  =0,  /72=2,  and  a=4°. 
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Figure  6.3.  Streamfunction  contours  at  equally  spaced  time  intervals  during  one 
complete  oscillation  of  the  freestream  velocity  for  the  case  of  the  previous  figure. 
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Figure  6.4.  Constant  pressure  contours  at  equally  spaced  time  intervals  during  one 
complete  oscillation  of  the  freestream  velocity  for  the  case  of  the  previous  figure. 
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Figure  6.5.  Time  series  of  system  response  for  an  elastic  membrane  wing  subjected  to 
harmonic  freestream  oscillation  is  shown  in  (a).  Surface  pressures  and  membrane 
profiles  at  several  times  during  one  complete  oscillation  of  the  freestream  velocity  are 
shown  in  (b)  and  (c),  respectively.  The  nondimensional  parameters  are  e=0,  Rn=4xl03, 
St=l.5,  III  =7.9,  U2  =0,  and  a=4°. 
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Figure  6.6.  Streamfunction  contours  at  equally  spaced  time  intervals  during  one 
complete  oscillation  of  the  freestream  velocity  for  the  case  of  the  previous  figure. 
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Figure  6.7.  Constant  pressure  contours  at  equally  spaced  time  intervals  during  one 
complete  oscillation  of  the  freestream  velocity  for  the  case  of  the  previous  figure. 
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The  stream  function  contours  for  this  case  are  shown  in  Fig.  6.6  and  the  associated 
pressure  contours  are  shown  Fig.  6.7.  It  may  be  seen  by  comparing  Fig.  6.3  and  Fig.  6.6  that 
although  the  two  cases  have  nearly  identical  steady  state  behavior  (recall  the  solution  at  t=0 
is  the  steady  state  solution  in  both  figures)  the  unsteady  response  of  the  two  cases  is 
dramatically  different.  The  differences  between  the  two  cases  may  be  broadly  attributed  to 
the  geometric  nonlinearity  inherit  in  the  elastic  problem. 

6.3  Inextensible  Membrane  Case 

The  final  class  of  problems  to  be  investigated  is  the  case  of  an  elastic  wing  with  the 
material  stiffness  tending  to  infinity.  The  controlling  dimensionless  parameter  8,  was  chosen 
to  coincide  with  a  portion  of  the  experimental  data  reported  by  Newman  and  Low  (1984). 
The  value  of  the  aeroelastic  parameter  IJj ,  was  chosen  large  enough  so  that  the  membrane 
may  be  considered  essentially  inextensible.  Figure  6.8  (a)  shows  the  time  series  of  the 
aerodynamic  and  membrane  properties  for  the  inextensible  wing  at  0°  angle  of  attack.  It  may 
be  seen  from  the  figure  that  the  X2  coordinate  of  the  membrane  at  the  midchord  point  is 
largely  unaffected  by  variation  in  the  freestream  velocity.  However  the  membrane  profile 
does  change  shape  in  response  to  the  unsteady  flow  field.  This  shape  change  is  shown  is  Fig. 
6.8  (c)  at  two  time  instants  during  the  second  complete  cycle  of  the  freestream  velocity.  For 
this  class  of  problems  the  membrane  profile  is  responding  primarily  to  the  periodic 
appearance  and  collapse  of  regions  of  recirculating  flow  on  both  the  upper  and  lower 
surfaces  of  the  membrane.  These  regions  of  recirculating  flow  may  be  seen  in  Fig.  6.9.  It  may 
be  seen  that  as  the  freestream  begins  to  decelerate  a  region  of  separated  flow  appears  on  the 
lower  surface  of  the  membrane  and  a  pair  of  counterrotating  recirculating  zones  appear  on 
the  upper  surface  of  the  membrane  near  the  trailing  edge.  As  the  freestream  velocity  begins 
to  accelerate  these  recirculating  zones  quickly  collapse  and  the  flow  once  again  becomes 
fully  attached.  As  was  the  case  with  the  constant  tension  wing  the  aeroelastic  response  shown 
in  Fig.  6.8  may  not  be  characterized  as  simple  harmonic  response  at  the  freestream  forcing 
frequency.  The  associated  pressure  contours  are  shown  in  Fig.  6.10. 
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Figure  6.8.  Time  series  of  system  response  for  an  inextensible  membrane  wing 
subjected  to  harmonic  freestream  oscillation  is  shown  in  (a).  Surface  pressures  and 
membrane  profiles  at  several  times  during  one  complete  oscillation  of  the  freestream 
velocity  are  shown  in  (b)  and  (c),  respectively.  The  nondimensional  parameters  are 
e=.017,  Rn=4xl03,  St=1.5, 77;  =46, 772=0,  and  a=0°. 
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Figure  6.9.  Streamfunction  contours  at  equally  spaced  time  intervals  during  one 
complete  oscillation  of  the  freestream  velocity  for  the  case  of  the  previous  figure. 
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Figure  6.10.  Constant  pressure  contours  at  equally  spaced  time  intervals  during  one 
complete  oscillation  of  the  freestream  velocity  for  the  case  of  the  previous  figure. 


CHAPTER  7 
A  THREE-DIMENSIONAL  MEMBRANE  WING  MODEL 

7.1  The  Elastic  Membrane  Problem  in  Three-Dimensions 

In  this  section  the  principle  of  virtual  work  is  used  to  formulate  the  finite  element 
matrices  for  a  bilinear  quadrilateral,  plane  stress  membrane  element.  The  total  Lagrangian 
formulation  is  used  in  deriving  the  element  which  may  undergo  large  displacements  and 
rotations  but  is  assumed  to  undergo  only  small  strains.  The  formulation  and  nomenclature 
used  follows  the  work  of  Bathe  (1984).  Index  notation  with  the  summation  convention  is 
used  throughout  and  a  comma  implies  partial  differentiation  with  respect  to  the  indicated 
independent  variable. 


configuration  at  time  t  +  At 


configuration  at  time  t 


Figure  7.1.  Deformation  trajectory  of  material  volume  °V 

Figure  7. 1  shows  a  material  volume  undergoing  large  deformation  from  a  natural 
unstrained  state  at  time  t=0  through  an  intermediate  state  at  time  t  to  a  final  deformed 
configuration  at  time  t+At.  The  time-like  coordinate  t  is  used  here  only  as  a  convenient 
parameter  for  describing  the  loading  and  deformation  state  and  does  not  imply  the  existence 
of  a  dynamic  process.  The  statement  of  equilibrium  for  the  continuum  at  configuration  t+A  t 
is  given  by  the  principle  of  virtual  work  as 
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l+AtSijd   {l+AtEij)odv  =   ,+A,R  (?1) 

J  oy 

where  Sy  is  the  2nd  Piola-Kirkoff  stress  tensor,  e,y  is  the  Green's  strain  tensor,  6  is  an 
admissible  variation  and  t+AtR  is  the  virtual  work  performed  by  external  forces  on  the 
material  volume  at  time  t+A  t. 

Consider  next  the  incremental  decomposition  of  the  stresses,  strains,  displacements, 
and  material  coordinates  given  by 

t+AtF  ,   _  tF      ,    „  (7.2) 

tij  cij    T    cy 

<^<S,   =  '5,  +  S,  (7-3) 

'+zl'M/   =  \.  +  ut  (7.4) 

'^'x,.   =  '*,-  +  ut  (7-5) 

where  £,y  ,  5y- ,  and  u\  are  increments  in  strain,  stress,  and  displacement,  respectively.  Using 
these  decompositions  in  Eq.  (7.1)  along  with  the  following  approximations 

Sij  ~  Cijrs  en  (7-6) 

d  By  ~  d  etj  (7-7) 

gives  the  following  linearized  incremental  equation  of  equilibrium 

'Sij  6  etj   °dV    (7.8) 


Cijrsersdeij°dV+   [      'SlJdrllJ°dV  =  l+*'R  - 

J   oy 


oy 


where  Qjrs  is  the  material  property  tensor,  and  etj  and  rjij  are  the  linear  and  nonlinear 
components,  respectively,  of  the  incremental  strain  tensor  given  by 

eH  =  2(  "<V  +  UJ*  +  '"*•<  UkJ  +  "*•«'  '"^'  }  (7'9) 
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Having  stated  the  general  incremental  equilibrium  condition  in  Eq.  (7.8)  consider 
the  development  of  the   finite  element  matrix  expressions   for  a  plane   stress, 
isoparametric,  bilinear  quadrilateral  element  shown  in  Fig.  7.2. 


now 


*2 


X3 


xi 


global  system 


local  system 


Figure  7.2.  Plane  stress  quadrilateral  element  with  local  and  global  coordinate  systems 

The   element   coordinates,    displacements,    and    displacement    increments    are 
interpolated  from  their  nodal  values  in  the  element  coordinate  system  as  follows 


4   iT 


=  [H][X\    X\    X\    X\.  .  .X\\ 


(7.11) 


=  [//][    'U\    'U\    'U\    'U]    .  .  .'U\]T 


(7.12) 


where 


=  [H][U\     U\    U\    U\ 


..u\\T 


[H]    = 


h\    °    °  h2  0  0  h3  0  0  hA  0  o' 
0   hx    0  0/i2°  0/j3  °  0^4° 
0    0    ft,  0  0  h2  0  0  h3  0  0  h\ 


(7.13) 


(7.14) 


with  the  isoparametric  interpolation  functions  hj,  hi,  /ij,  and  I14  given 
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ht   =  \  (  1  -  r)(  1  -  s)  (7.15) 


i 
4 


/*2    =  |  (  1   +   r  )  (  1   -  5  )  (7.16) 


h3    =  I  (  I    +   r)  (  I    +   s)  (7.17) 

A4=i(l-r)(l+j)  (7.18) 

where  r  and  5  are  natural  coordinates  defined  on  the  bi-unit  square.  Substituting  these 
coordinates  and  displacements  interpolations  into  Eq.  (7.8)  and  invoking  the  principle  of 
virtual  work  for  each  nodal  displacement  in  turn  produces  the  following  set  of  discrete 
equilibrium  equations  in  index  free  form 

(  [  lKL  ]   +   [  *KNL  ]  )  [  U  ]  «    =   [  t+AtR  ]  ('"1)    -   [  t+AtF  ]  G-1)  (7-19) 

where  the  superscript  indicates  the  iteration  level  of  the  solution  procedure.  The  element 
stiffness,  force,  and  load  matrices  appearing  in  Eq.  (7.19)  are  defined  as  follows 


'KL]=    (     ['i 


['KL]=  [!BL]'  [  C]  ['BL]°dV  (7.20) 

"V 

\-'KNL]=\      [!BNL]TVS}[!BNL}°dV  (7.21) 

[t+AtF]=    '      [<+*'BL]T[!+A,S]°dV  (7.22) 

J  oy 
l  +  At  n    i      _       I         r     tj   -\T   r  t  +  At  r  i  o 


I- 

J  °S 


[,+n,R]    =    ;      [  H]'  [t+Atf]°dS  (7.23) 

°s 

where  the  explicit  expressions  for  both  the  linear  and  nonlinear  strain  displacement  matrices, 
[  Bi  ]  and  [  B^i  ],  respectively,  may  be  found  in  Bathe  (1984). 

At  this  point  a  departure  is  made  from  the  general  three-dimensional  continuum 
formulation  by  assuming  a  state  of  plane  stress  to  exist  in  the  element  and  employing  a  plane 
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stress  consitutive  law.  These  assumptions  lead  to  the  following  choice  for  the  material 
property  matrix 


[  C]    = 


1  -  v2 


1 

V 

0 

V 

1 

0 

0 

0 

1 

— 

V 

2 

,where  v  is  Poisson's  ratio,  and  for  the  stress  state  matrices 


(7.24) 


Vs  ]  = 


1 S      "\        0       0 

^  11      ^12     u         u 

,r      re        0        0 

15 12      ^22    tr        !r 

o     o     *» 


0 
0 

0 


0  '5i2  '5: 
0  0  0 
0       0       0 


0 

0    ' 

0 

0 

12     0 

0 

22  tv 

0 
0 

°12 

'sn 

°22 

(7.25) 


and 


t+At 


[,+"'S  ]    = 


t+At 


t+At 


t+At 


'22 
>12 


(7.26) 


With  these  choices  for  the  constitutive  and  stress  state  matrices  the  linear  strain 
transformation  matrix  becomes 


where 


[<BL]    =   ['BL0]  +  [lBn  ] 


(7.27) 


I'Brn]    = 


dh , 

dx. 


0  0 


dh2 

ax7 


dh, 

0^00 


dhx  dhx 


0 


dh* 


dx2    ^x  1    ^x2 


.  0 
.  0 
.  0 


(7.28) 


and 


[  'B,A  = 


with 


/ 


/, 


"ax"; 

dh, 
2~dx, 


•11 


dh, 
dx~ 


+  I 


dh. 
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dh. 


dh, 
22  3x7 

+  L 


dh. 


dh, 

'si  a^ 

a/2, 


a/z, 


12  ax,     '21ax,     '22ax,     '31  ax 


+  / 


dh. 


BhA 

'31  ax, 
,  dh4 
'32  aZ 


a/i, 


32  ax,  •  *3i3jc„ 


+  / 


a/?, 


32ax, 


(7.29 


a^ 


H«ri* 


/t>-   =   hrr    '£/•        U=  1,2,3    and    fc  =  1,2,3,4 


(7.30) 


Similarly,  the  nonlinear  strain  transformation  matrix  becomes 


[  'BNL   ] 


a/z, 
ax, 

d/t, 

ax0 


o 
o 
o 
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0      0 
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dh, 
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a/i3 
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dh, 
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(7.31) 


Since  the  elastic  membrane  model  described  above  will  be  used  in  conjunction  with 
an  inviscid  flow  model  the  surface  traction  vector  is  taken  to  be  normal  to  the  element  and 
equal  to  the  pressure  difference,  Ap,  across  the  membrane.  This  pressure  difference  is 
assumed  to  be  uniform  over  each  element.  Consequently,  the  surface  traction  vector  is  given 
by 
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t+At 


V™f]     = 


[01 

0        ] 

0 

— 

0 

?p. 

p 

-p\ 

(7.32) 


where  p+  and/?-  are  the  aerodynamic  surface  pressures  on  the  upper  and  lower  membrane 
surfaces  respectively. 

The  transformation  relating  the  derivatives  in  physical  coordinates  to  the  derivatives 
in  isoparametric  coordinates  is  given  by  the  chain  rule.  This  transformation  is 


[  d  1 

dXy 

d 

dx2 

■&) 

dx2 
17 


dx2 
~dr 


dx,       dx. 


ds 


dr 


dr 
ds 


(7.33) 


where 


J  = 


dX]  dx2       3jC]  dXf 


(7.34) 


/■ 


dr   ds         ds   dr 

The  volume  integrals  appearing  in  Eqs.  (7.20)  through  (7.23)  are  evaluated  using 
3x3  Gauss  quadrature  on  the  bi-unit  square  according  to 

+  i    +i 
(  )°dV  =  t  J  (  )  dr  ds  (7.35) 

nv  -i-i 

where  t  is  the  thickness  of  the  element. 

Since  the  element  matrices  were  evaluated  in  a  local  element  coordinate  system  the 
element  matrices  must  be  transformed  to  the  global  coordinate  system  prior  to  assembly. 
These  transformations  are  accomplished  using  the  conventional  1st  and  2nd  order  tensor 
transformations  given  by 

,7  r    tr  i    r    <r  i  (7-36) 


[K]=[T]'[K][T] 


[F]    =   [  T]T  [  F] 


(7.37) 


[R]    =  [T]T[R] 


(7.38)  8) 
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where  [  T  ]  is  the  matrix  of  direction  cosines  relating  the  global  system  to  the  local  system 
(Malvern  1969).  After  transformation  the  system  of  global  equilibrium  equations  is 
assembled  using  the  direct  stiffness  method  Bathe  (1984).  The  nodal  displacement 
increment  residual  may  be  defined  as  follows 

|  — —  (7.39) 

alldof 

where  c  is  the  membrane  chord  and  the  sum  is  over  all  unconstrained  degrees  of  freedom. 

It  should  be  noted  that  the  linearizing  approximations  given  by  Eqs.  (7.6)  and  (7.7) 
require  that  an  iterative  procedure  be  used  to  satisfy  the  exact  equilibrium  condition  given 
in  Eq.  (7.1).  The  procedure  described  above  for  incrementally  increasing  the  aerodynamic 
pressure  on  the  membrane  and  iterating  until  a  converged  solution  is  achieved  provides 
considerable  flexibility  in  the  choice  of  a  solution  strategy  for  aeroelastic  problems.  The 
particular  strategy  adopted  for  the  membrane  wing  model  will  be  discussed  in  the  section 
on  computational  results. 

A  very  appealing  feature  of  the  finite  element  formulation  described  above  is  that  the 
assembled  stiffness  matrix  is  symmetric  and  narrowly  banded.  This  feature  is  largely 
responsible  for  the  usefulness  of  the  method  in  solving  problems  where  the  number  of 
unknowns  is  large.  This  is  in  direct  contrast  to  the  situation  encountered  in  the  aerodynamic 
problem  where  the  assembled  influence  coefficient  matrix  is  neither  symmetric  nor  banded. 
Consequently  as  the  number  of  aeroelastic  elements  increases  the  aerodynamic  problem 


requires  an  increasingly  larger  portion  of  the  total  solution  time. 

7.2  The  Aerodynamic  Problem  in  Three-Dimensions 

In  this  section  a  method  is  presented  for  determining  the  forces  and  moments  as  well 
as  the  pressure  distribution  resulting  from  the  inviscid  flow  over  thin  wings  of  finite  span. 
The  fundamental  assumption  concerning  the  flow  field  around  the  wing  is  that  it  is 
irrotational,  and  consequently  the  velocity  field  may  be  derived  from  a  scalar  potential. 
Implicit  in  this  assumption  is  that  the  flow  is  inviscid,  incompressible,  and  free  of  vorticity 
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far  upstream.  The  use  of  vortex  filament  singularities  to  model  the  lifting  potential  flow 
around  thin  wings  has  its  historical  origin  and  justification  in  work  by  James  (1972).  A 
description  of  the  method  in  its  modern  form  for  wings  of  finite  span  may  be  found  in  Katz 
and  Plotkin  (1991). 
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Figure  7.3.  Discretization  of  a  elliptical  planform  thin  wing  into  a  number  of 
quadrilateral  vortex  filament  panels. 


Figure  7.3  shows  a  thin  wing  that  has  been  discretized  into  a  number  of  quadilateral 
panels.  A  typical  panel  is  composed  of  a  horseshoe  shaped  vortex  filament  and  a  control  point 
where  the  flow  tangency  condition  is  enforced.  The  figure  also  shows  the  vortex  filament 
to  extend  downstream  to  infinity  as  required  by  the  Kelvin-Helmholtz  vorticity  theorem 
(Katz  and  Plotkin  1991).  The  existence  of  this  trailing  vortex  wake  is  the  unique  feature  of 
a  finite  span  wing  that  distinguishes  it  from  a  two-dimensional  airfoil. 
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By  modelling  the  wing  as  an  assemblage  of  vortex  singularity  segments  of  unknown 
strength  and  enforcing  the  zero  normal  relative  velocity  condition  at  each  control  point  the 
following  set  of  linear  algebraic  equations  may  be  formed  and  solved  for  the  vortex  strengths 

r 

t+At[  [A][r]=[-v.  •«]]  (7.40) 

where  [  A  ]  is  the  matrix  of  vortex  influence  coefficients,  v  «,  is  the  freestream  velocity  vector, 
and  n  is  the  panel  unit  normal  vector.  If  the  vortex  filaments  and  control  points  are  located 
at  the  local  quarter  and  three-quarter  chord  points  of  the  panel,  respectively,  the  Kutta 
condition  is  implicitly  satisfied  and  no  additional  boundary  condition  is  needed  at  the  trailing 
edge  (James  1972,  Katz  and  Plotkin  1991).  For  the  sake  of  brevity  the  superscript  in  Eq. 
(7.40)  indicating  the  deformation  state  will  be  omitted  in  the  following  development  of  the 
aerodynamic  problem  since  it  is  clear  that  the  aerodynamic  pressures  are  to  be  evaluated  in 
the  final  equilibrium  state  at  time  t+A  t. 

Once  the  circulation  strength  for  each  panel  has  been  determined  from  Eq.  (7.40)  the 
force  vector,  F  acting  on  the  bound  segment  of  each  vortex  filament  is  given  by  the 
generalized  Kutta-Joukowski  theorem  (Katz  and  Plotkin  1991)  as 

F   =  Ql  (  v   x  r  )  (7.41) 

where  q  is  the  fluid  density,  /  is  the  length  of  the  bound  segment  of  the  vortex  filament,  v  is 
the  velocity  vector  at  the  control  point,  and  f  is  the  circulation  vector  per  unit  length  of  the 
bound  vortex  filament.  Having  determined  the  net  force  acting  on  the  bound  vortex  filament 
using  Eq.  (7.41),  the  net  pressure  difference  acting  on  the  aeroelastic  element  is  then  simply 
given  by 

Ap   =  |  (7.42) 

where  F  is  the  magnitude  of  F  and  a  is  the  panel  area. 
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Figure  7.4.  Resolution  of  far  field  panel  forces  into  lift  and  drag  components 

The  force  acting  on  each  panel  is  again  determined  by  using  the  generalized 
Kutta-Joukowski  theorem  except  now,  in  order  to  account  correctly  for  the  downwash  and 
induced  drag  on  the  wing,  the  velocity  that  is  used  v*,  is  the  sum  of  the  freestream  velocity 
and  the  wake  induced  velocity  (Katz  and  Plotkin  1991).  Figure  7.4  illustrates  the  far  field 
panel  force  resolved  into  lift  and  drag  components  and  are  seen  to  be  given  by 


L   =  F3Cosa  —  Fjsina 


(7.43) 


D   =  f^sina  +  Fjcosa 


where  a  is  the  angle  of  attack  and  the  far  field  panel  force  is  given  by 


(7.44) 


F*    =  gl  (  v*    x   r  )  (7.45) 

The  total  lift  and  drag  on  the  wing  is  obtained  by  summing  up  the  contributions  from  each 
bound  vortex  filament  over  the  entire  wing. 

The  aerodynamic  analysis  procedure  described  above  offers  a  simple  method  for 
approximating  the  aerodynamic  characteristics  of  thin  wings  of  finite  span.  Unfortunately 
the  potential  flow  model  is  a  very  restrictive  fluid  dynamic  model.  Many  important  features 
of  real  fluid  behavior  such  as  boundary  layer  growth  and  separation  are  unapproachable  from 
a  potential  flow  viewpoint.  However  the  potential  flow  model  does  serve  as  a  reliable  tool 
for  developing  a  computational  procedure  for  the  aeroelastic  analysis  of  membrane  wings 
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of  finite  span.  The  synthesis  of  the  plane  stress  quadrilateral  element  and  the  quadrilateral 
vortex  panel  to  form  a  quadrilateral  membrane  wing  element  is  described  in  the  next  section. 

7.3  The  Aeroelastic  Problem 

The  development  of  a  computational  procedure  for  the  aeroelastic  analysis  of 
membrane  wings  may  be  divided  into  two  parts  -  the  development  of  a  viable  membrane 
wing  element  and  the  development  of  a  solution  strategy  for  the  coupled  aeroelastic 
boundary  value  problem.  In  developing  such  a  procedure  there  are  many  options  available 
in  terms  of  both  the  aeroelastic  element  technology  and  the  solution  strategy  adopted  for  the 
solution  of  the  coupled  boundary  value  problem.  Of  course,  the  problem  formulation  and 
discretizations  that  were  described  in  the  previous  sections  were  chosen  such  that  a  hybrid 
aeroelastic  element  could  be  easily  formed  by  combining  the  plane  stress  quadrilateral 
element  of  section  7. 1 .  with  the  quadrilateral  vortex  filament  panel  of  section  7.2. 

Recalling  the  load  stepping  procedure  described  in  section  7. 1  for  solving  the  elastic 
problem,  it  is  natural  to  extend  this  methodology  to  a  velocity  stepping  procedure  for  the 
solution  of  the  aeroelastic  problem.  Furthermore,  since  a  number  of  iterations  are  required 
at  each  load  step  in  the  elastic  problem  in  order  to  satisfy  equilibrium,  it  is  again  natural  to 
evaluate  the  aerodynamic  pressures  based  on  the  updated  wing  configuration  during  each 
iteration.  This  procedure  guarantees  that  the  wing  configuration  and  the  aerodynamic 
loading  are  in  equilibrium  (within  some  prespecified  tolerance)  at  the  end  of  each  velocity 
step.  However,  iterating  at  each  load  step  in  order  to  achieve  aeroelastic  equilibrium  is  quite 
expensive  since  the  computation  of  the  aerodynamic  surface  pressures  involves  the  solution 
of  a  large  nonbanded  system  of  equations.  Consequently,  the  CPU  time  required  to  solve  the 
aerodynamic  problem  increases  much  more  rapidly  with  problem  size  than  the  CPU  time 
required  to  solve  the  elastic  problem.  This  situation  is  a  direct  result  of  the  handedness  of 
the  stiffness  matrix  and  the  fullness  of  the  aerodynamic  coefficient  matrix. 

The  iterative  aeroelastic  procedure  imbedded  within  a  velocity  stepping  procedure 
as  described  above  is  a  very  general  strategy  that,  while  guaranteeing  aeroelastic  equilibrium 
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at  each  velocity  step,  may  not  be  the  most  efficient  or  reliable  strategy  for  solving  certain 
membrane  wing  problems.  Indeed  the  strategy  may  fail  to  converge  in  some  situations. 
Consequently  some  tailoring  of  the  strategy  may  be  necessary  in  some  cases.  Some  of  these 
situations  will  be  described  in  the  next  section  on  computational  results. 

7.4  Three-Dimensional  Test  Cases 

Before  applying  the  procedure  described  in  the  previous  sections  to  the  analysis  of 
membrane  wings  two  limiting  cases  of  the  aeroelastic  boundary  value  problem  are 
investigated  and  the  results  of  the  computational  procedure  are  compared  with  well  known 
analytic  solutions.  The  two  limiting  cases  presented  in  this  section  are  the  deformation  of 
a  square,  edge  constrained,  elastic  membrane  subjected  to  uniform  pressure  load  and  the 
aerodynamic  characteristics  of  a  rigid  ,  zero  thickness  elliptical  planform  wing  of  moderate 
aspect  ratio. 
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Figure  7.5.  Comparison  of  computed  and  analytic  midpoint  displacement,  (b),  of  a 
square,  initially  flat  membrane  subjected  to  uniform  pressure  load.  Computed  solution 
used  256  square  elements  as  shown  in  (a).  Poisson's  ratio  is  0.3. 

Figure  7.5  (a)  shows  an  initially  flat,  square,  edge  constrained  membrane  that  has 
been  discretized  into  256  square  plane  stress  elements.  Figure  7.5  (b)  shows  the  computed 
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transverse  midpoint  deflection  of  the  membrane  subjected  to  a  uniform  pressure  load. 
Results  are  shown  in  the  figure  for  several  different  values  of  the  dimensionless  stiffness 
parameter  III,  where  uniform  pressure  has  been  substituted  for  the  stagnation  pressure  in  the 
definition  of  rij.  The  analytic  result  for  this  problem  is  shown  by  the  solid  line  in  the  figure 
and  is  taken  from  Seide  (1977).  The  agreement  between  the  computed  solution  and  the 
analytic  solution  is  seen  to  be  quite  good  even  with  such  a  a  small  number  of  elements.  For 
completeness,  all  three  components  of  the  displacment  field,  are  shown  in  Fig.  7.6  where  the 
displacement  profiles  are  taken  along  lines  where  the  X2  coordinate  is  constant. 

Since  the  membrane  is  initially  flat  and  free  of  any  pretension,  the  tangent  stiffness 
matrix  is  singular  for  the  problem  of  an  initially  flat  membrane  without  pretension,  and  the 
load  stepping  procedure  fails  for  the  first  load  step.  Consequently,  the  solution  strategy 
adopted  for  the  test  case  of  Figs.  7.5  and  7.6  was  to  pretension  the  membrane  for  the  first 
load  step  in  order  to  remove  the  singularity  in  the  tangent  stiffness  matrix,  and  subsequently 
remove  the  pretension  before  assemblying  the  tangent  stiffness  at  the  next  load  step. 

The  solution  shown  in  the  Fig.  7.6  was  abtained  by  taking  5  load  steps  to  develop  the 
final  pressure  load  and  allowing  iteration  only  at  the  last  load  step  when  the  pressure  had 
reached  its  final  value.  The  convergence  path  for  this  problem  may  be  seen  in  Fig.  7.6  (d) 
It  may  be  seen  that,  once  the  iterative  procedure  begins  after  load  step  5,  the  covergence  is 
very  rapid.  However,  if  fewer  load  steps  were  taken  before  iteration  was  allowed,  the 
convergence  would  have  been  much  slower  since  the  tangent  stiffness  matrix,  at  load  step 
2  or  3,  would  have  been  a  much  less  accurate  estimate  of  the  stiffness  of  the  membrane  than 
at  load  step  5.  Typically,  between  5  and  10  explicit  load  steps  are  adequate  to  insure  a  stable 
and  rapid  convergence  path  for  uniform  pressure  loading.  However  for  an  aeroelastic 
problem,  more  explicit  load  steps  may  be  required  prior  to  iteration  since  the  net  pressure 
acting  on  the  membrane  will  then  depend  on  the  membrane  configuration. 
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Figure  7.6.  Computed  displacements  of  a  square,  edge  restrained,  initially  flat 
membrane  subjected  to  uniform  pressure  load  with  77/  =  2.5.  Computed  solution  used 
256  square  elements.  Poisson's  ratio  is  0.3. 


Having  investigated  the  accuracy  of  the  elastic  formulation  with  a  uniform  pressure 
load  in  the  previous  section  consider  now  the  aerodynamic  properties  of  a  thin,  rigid, 
cambered,  elliptic  planform  wing  in  uniform  flow.  Figure  7.7  shows  the  computed 
chordwise  loading  at  several  spanwise  locations  for  an  elliptical  wing  of  aspect  ratio  10  with 
the  xj  coordinate  defined  by  a  NACA  6-series  meanline  (Abbott  and  von  Doenhoff  1959). 
The  angle  of  attack  was  adjusted  until  the  difference  between  the  geometric  angle  and  the 
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induced  angle  was  equal  to  the  ideal  angle  of  attack  for  the  chosen  meanline.  In  Fig.  7.7  the 
computed  chordwise  loading  at  several  spanwise  locations  along  the  wing  is  shown  along 
with  the  two-dimensional  analytic  result  from  thin  wing  theory  for  the  same  NACA 
meanline  at  the  ideal  angle  of  attack  (Abbott  and  von  Doenhoff  1959).  The  computed  loading 
on  the  midspan  cross  section  may  be  seen  to  be  in  good  agreement  with  the  analytic  result 
except  very  near  the  leading  edge  where  the  potential  flow  model  is  singular. 
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Figure  7.7.  Computed  pressure  coefficients  at  several  spanwise  locations  for  a  rigid, 
ellipitical  planform  wing  of  aspect  ratio  10  operating  at  the  ideal  angle  of  attack 
(a=2.58°).  The  profile  of  the  wing  is  defined  by  a  NACA  6-series  meanline  (a=.50, 
Cli=-50).  The  computed  solution  used  15  spanwise,  and  30  chordwise  panels. 

The  accuracy  of  the  lifting  surface  model  may  be  further  tested  by  computing  the  lift 
curve  slope,  C^a,  for  an  elliptic  planform  wing  of  varying  aspect  ratio.  A  comparison  is 
shown  in  Fig.  7.8  for  three  different  aspect  ratios  between  3  and  10.  The  analytic  result  is 
taken  from  Van  Dyke  (1975).  The  agreement  with  the  analytic  result  may  be  seen  to  be  quite 
good. 
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Figure  7.8.  Comparison  of  computed  and  analytic  lift  curve  slopes  for  an  elliptical 
planform  wing  of  varying  aspect  ratio.  Computed  solution  used  24  spanwise,  and  12 
chordwise  panels  for  a  total  of  288  quadrilateral  elements. 


7.5  Application  to  an  Elliptical  Planform  Elastic  Wing 

Having  established  the  accuracy  of  the  elastic  and  aerodynamic  computational 
procedures  independently  in  the  previous  section  consider  next  the  application  of  the 
procedure  to  a  flexible  membrane  wing.  Figure  7.9  shows  the  computed  equilibrium 
chordwise  loading  and  membrane  geometry  for  an  initially  flat,  pretensioned,  edge 
constrained,  elliptic  planform  wing  of  aspect  ratio  4.  The  controlling  aeroelastic  parameter 
is  772  and  is  taken  to  be  equal  to  2  and  the  angle  of  attack  is  5°.  For  comparison,  the  chordwise 
loading  for  rigid  wing  of  the  same  geometry  and  operating  under  the  same  conditions  is  also 
shown  in  Fig  7.9  (b).  As  may  be  seen  in  the  figure  the  effect  of  stretch  in  the  membrane  is 
to  move  the  center  of  pressure  toward  the  trailing  edge  of  the  wing  and  to  increase  the  lift. 
The  velocity  stepping  strategy  used  for  this  problem  was  to  take  10  steps  to  develop  the  final 
velocity  and  to  allow  iteration  only  at  the  last  load  step.  The  convergence  path  shown  in  Fig. 
7.9  (d)  is  at  the  10th  load  step. 
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Figure  7.9.  Equilibrium  chordwise  aerodynamic  loading  for  (a),  an  initially  flat, 
pretensioned,  ellipitical  planform,  constant  tension  membrane  wing  of  aspect  ratio  4,  77 
=  2, 77;  =  0,  and  a  =  5°.  and  (b),  a  similar  rigid  wing  at  the  same  angle  of  attack,  (c), 
equilibrium  membrane  profiles,  (d),  residual  of  incremental  nodal  displacements. 
Poisson's  ratio  is  0.3.  Computed  solutions  used  15  spanwise,  and  30  chordwise  panels. 

The  effect  of  membrane  deformation  on  wing  lift  can  be  seen  in  Fig.  7.10  where  the 
lift  curves  for  both  a  constant  tension  wing  and  an  elastic  wing  are  shown  along  with  the  lift 
curve  for  a  rigid  wing  of  the  the  same  aspect  ratio.  The  effect  of  membrane  stretch  on  the 
constant  tension  wing  is  simply  to  increase  the  lift  curve  slope  by  an  amount  inversely 
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proportional  to  77?,  otherwise  the  variation  of  lift  with  angle  of  attack  is  identical  to  a  rigid 
wing.  The  behavior  of  the  lift  curve  for  the  elastic  wing  is  somewhat  different  due  the 
intrinsic  nonlinearity  in  the  deformation  of  the  membrane.  At  very  small  incidence  angles 
the  elastic  wing  develops  lift  very  rapidly  since  the  membrane  cannot  initially  resist 
deformation.  However  once  the  wing  inflates  and  geometrically  stiffens  there  is  very  little 
further  deformation  and  the  lift  curve  slope  is  nearly  identical  to  a  rigid  wing  of  the  same 
aspect  ratio. 
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Figure  7.10.  Variation  of  lift  coefficient  with  angle  of  attack  for  initially  flat,  ellipitical 
planform,  elastic  and  constant  tension  membrane  wings  of  aspect  ratio  4.  Poisson's  ratio 
is  0.3.  Computed  solutions  used  16  panels  spanwise  and  chordwise. 


Finally,  the  sensitity  of  the  aeroelastic  solution  to  the  discretization  scheme  is  shown 
in  Table  7. 1 .  As  may  be  seen  in  the  table  the  computed  lift  is  essentially  the  same  for  all  three 
discretizations  tabulated.  The  computed  drag  is  somewhat  more  sensitive  to  the  refinement 
of  the  grid  since  it  is  proportional  to  square  of  the  lift  coefficient.  However  the  relative 
insensitivity  of  results  to  the  discretization  is  largely  an  artifact  of  the  edge  constrained 
boundary  condition  used  in  the  computations.  For  other  edge  boundary  conditions,  e.g..  the 
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free  trailing  edge  of  a  marine  sail,  the  computed  solution  would  undoubtedly  demonstrate 
a  greater  sensitivity  to  grid  refinement. 


Table  7.1.  Effect  of  grid  refinement  on  the  computed  lift  and  drag  for  an  initially  flat, 
ellipitical  planform,  elastic  membrane  wing  of  aspect  ratio  4.  TIj  =  10,  IJ2  =  0,  and  a  = 

5°.  Poisson's  ratio  is  0.3. 
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CHAPTER  8 
SUMMARY  AND  CONCLUSIONS 

In  the  present  work  a  numerical  model  simulating  the  aeroelastic  characteristics  of 
a  flexible  two-dimensional  membrane  wing  has  been  presented.  The  use  of  the 
Navier-Stokes  equations  as  the  fluid  dynamic  model  in  the  present  model  is  a  substantial 
departure  from  previous  work  on  sail  aerodynamics  which  has,  almost  universally,  adopted 
a  potential  based  description  of  the  flow  field. 

The  two-dimensional  aeroelastic  boundary  value  problem  was  nondimensionalized 
and  a  set  of  six  basic  dimensionless  parameters  was  derived  which  govern  the  solution  of  the 
problem.  An  additional  parameter,  the  frequency  ratio,  was  proposed  as  a  meaningful 
parameter  for  characterizing  the  harmonically  driven  unsteady  aeroelastic  problem. 

A  numerical  procedure  was  then  developed  for  the  solution  of  the  coupled  aeroelastic 
problem  and  was  shown  to  yield  results  that  are  in  agreement  with  available  analytic 
solutions  for  several  appropriate  limiting  cases.  The  numerical  procedure  was  also  shown 
to  satisfy  certain  identities  as  dictated  by  the  fundamental  fluid  dynamic  conservation  laws. 

The  role  of  viscosity  in  membrane  wing  aerodynamics  was  investigated  using  the 
numerical  model  for  both  steady  and  unsteady  flows.  These  investigations  were  facilitated 
by  distinguishing  three  classes  of  problems  which  are  associated  with  corresponding  limiting 
cases  of  the  dimensionless  parameter  set.  The  aerodynamic  characteristics  of  membrane 
wings  at  Reynolds  numbers  between  103  and  104  were  shown  to  be  substantially  different 
from  those  predicted  by  a  potential  based  membrane  wing  theory. 

The  role  of  viscosity  was  shown  to  be  preeminent  in  the  harmonically  forced 
unsteady  flow  about  a  membrane  wing.  In  this  case  the  influence  of  viscosity  is  enhanced 
since  the  acceleration  and  deceleration  of  the  freestream  velocity  strongly  influences  the 
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separation  and  reattachment  of  the  flow.  The  periodic  appearance  and  collapse  of 
recirculation  zones,  along  with  an  attendant  adjustment  in  the  membrane  configuration, 
results  in  an  aeroelastic  response  which  may  not  be  characterized  as  simple  harmonic 
response  at  the  freestream  forcing  frequency. 

Finally,  a  three-dimensional  potential  based  membrane  wing  element  was  derived 
and  shown  to  yield  results  that,  for  several  limiting  cases,  are  in  agreement  with  available 
analytic  solutions.  This  finite  span  aeroelastic  model,  although  it  makes  use  of  a  simplified 
fluid  dynamic  model,  nonetheless  serves  as  a  useful  first  step  in  developing  a  viscous  flow 
based  model  for  finite  span  membrane  wings. 

In  terms  of  future  work,  the  extension  of  the  present  two-dimensional  model  to 
include  the  effects  of  turbulence  in  the  flow  is  very  appealing.  This  extension  is  especially 
attractive  since  a  direct,  meaningful  comparison  of  the  numerical  predictions  could  then  be 
made  with  available  experimental  data.  Reasonable  agreement  between  the  measured  and 
predicted  aeroelastic  characteristics  of  membrane  wings  over  the  full  range  of  aeroelastic 
parameters  would  mark  a  turning  point  in  the  state  of  affairs  described  by  Marchaj  (1979) 
in  the  introduction. 
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